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INTRODUCTION 

A major problem in computational fluid dynamics is the generation of grids for 
realistic three-dimensional bodies such as complete aircraft configurations. The 
generation of a single grid that discretizes the entire flow region associated with 
a complex configuration is extremely difficult and, for some problems, it is imprac- 
tical. Often, an attempt to generate a single grid for a complicated flow region 
results in highly skewed grids which, in turn, result in inaccurate calculations; 
also, such an attempt usually requires an inordinate amount of the aerodynamicist ' s 
time to "tune" existing grid generation schemes to yield acceptable grids. A second 
factor that contributes to the complexity of grid generation is the necessity to 
cluster grid points in regions where the dependent variables and their gradients 
change rapidly (selective grid refinement). 

The problems mentioned above can be overcome, to a limited extent, by develop- 
ing very sophisticated grid generation schemes. However, an alternative approach is 
to divide the given region into simpler subregions such that each subregion (or 
zone) is a geometrically simple figure, e.g., every two-dimensional region can be 
divided into simple four-sided zones, and every three-dimensional region into six- 
sided zones. Grids can then be independently generated for each zone using existing 
grid-generation schemes. The "zonal" approach has the following advantages: 

1. The grid for any arbitrary region can be generated in a simple, straight- 
forward manner (since the zones are geometrically simple). 

2. Flow regions requiring grid refinement can be isolated in separate zones, 
and the required number of grid points can be introduced in these zones ( if neces- 
sary, this can be done adaptively as the solution progresses in time). 

3 . The zonal approach facilitates the use of different equation sets in dif- 
ferent zones; hence, simpler equation sets can be used in certain regions of the 
flow. This may result in a saving of computer time. 

4. The zonal approach also facilitates a block-processing scheme wherein only 
the data corresponding to certain regions of the flow field are required to reside 
in the main memory of the computer; the remaining data can be stored on disc or 
tape. Theoretically, the block processing technique permits the use of unlimited 
global grid sizes (computing speeds become a limiting factor in this case). 
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The subregions that result from the zoning process may be patched together 
(patched-grids) as in figure 1(a) or overlap one another (overlaid grids) as in 
figure 1(b). Both approaches have their relative advantages and disadvantages. 

Some disadvantages of overlaid grids are: a) A problem in n spatial dimensions 

requires interpolation in n dimensions (in order to transfer information from one 
grid onto another) while patched grids require only an (n - 1) dimensional interpo- 
lation (the details are given in a later section), b) Maintaining global conserva- 
tion seems to be more difficult with overlaid grids, c) The accuracy and conver- 
gence speed of the calculation seems to depend on the degree of overlap of the zones 
and the relative size of each zone, thus introducing a certain amount of undesirable 
empiricism in the formulation. However, the fact that zones can (and should) over- 
lap, may result in a certain amount of flexibility in generating grids in three 
dimensions. A second advantage of overlaid grids may lie in the ease with which 
grids can be moved relative to each other (for calculations involving bodies that 
move relative to each other). This set of lecture notes will deal primarily with 
patched-grids. 

The following examples serve to illustrate some of the advantages of the 
patched-grid approach. Figure 2 shows the flow region associated with a combination 
of three airfoils. The region is multiply connected and hence difficult to discre- 
tize with a single grid system. The division of the flow field into thirteen 
patches as in figure 2 results in simple four-sided regions which can be discretized 
easily. Figure 3(a) shows a six-finned projectile and figure 3(b) the grid used to 
discretize the flow region at the aft end of the projectile (ref. 1). Immediately 
after this axial station, an entirely new type of grid (fig. 3(c)) is required 
because of the sudden termination of the body. This transition can be achieved by 
treating the plane corresponding to the aft end of the projectile as a two- 
dimensional (surface) patch-boundary separating two three-dimensional patches. It 
should be noted that the grid in figure 3(c) is itself composed of four two- 
dimensional grids separated by one-dimensional (line) patch-boundaries. The grids 
in figure 3(c) alleviate the geometric singularity problem associated with a simple 
polar grid, and also facilitate the addition of grid points to the central portion 
of the flow region (to accurately model the wake). 

The division of a given region into patches introduces new boundaries in the 
calculation: the patch boundaries. Since the grid for each region is generated 

independently, the grid lines of two adjoining regions may align (continuous grids) 
or may not align (discontinuous grids) with each other. Even in the case of contin- 
uous grids, a sudden change in grid spacing or grid line orientation across the 
zonal boundary may give rise to discontinuities in the transformation metrics 
(metric-discontinuous grids). Figure 4 shows the different types of grids mentioned 
above . 

In order that information be transferred from one patch to another accurately, 
it is important to treat grid points on the patch boundaries with care. The nonlin- 
ear nature of the Euler and Navier-Stokes equations permits solutions with discon- 
tinuities, such as shocks. It is imperative that the finite-difference scheme used 
for the calculation be conservative so that these discontinuities, when captured, 
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assume the right strength and physical location. In a patched-grid calculation, it 
is important that the patch boundaries are also treated in a conservative manner so 
that the discontinuities can move freely across these boundaries. Results demon- 
strating the importance of maintaining conservation at patch boundaries can be found 
in references 2 and 3. 

Earlier work in the area of patched-grid calculations includes that of Cambier 
et al. (ref. 4). The authors analyze the patch-boundary problem for a system of 
hyperbolic equations and use the compatibility equations to develop a patch-boundary 
scheme. They present encouraging results for transonic channel flow. The use of 
the compatibility equations at the patch-boundary, however, results in a patch- 
boundary scheme that is not conservative; and hence, such a procedure is unsuitable 
for problems in which flow discontinuities move from one patch to another. 

Hessenius and Pulliam (ref. 2) present a conservative patch-boundary scheme for 
implicit integration methods like the Beam-Warming method (ref. 5).' Their results 
stress the need for a conservative patch-boundary scheme. However, the authors do 
not address the problem of discontinuous metrics or discontinuous grids across 
patch-boundaries. Rai, Chakravarthy and Hessenius (ref. 3) present results obtained 
on metric-discontinuous grids. The integration scheme used is the Osher upwind 
scheme. The proper choice of transformation metrics in the calculation of the 
fluxes makes the scheme fully conservative at both interior and patch-boundary 
points . 

A fully conservative patch-boundary scheme, that permits the movement of dis- 
continuities across patch boundaries with minimal distortion of the solution, is 
developed in reference 6. The scheme is designed for both discontinuous and metric- 
discontinuous grids and can be used in conjunction with first-order accurate, 
explicit integration schemes. The scheme is stable, accurate and is formulated in 
generalized coordinates. Results demonstrating these characteristics of the scheme 
are also presented in reference 6. The demonstration calculations include super- 
sonic flow over a cylinder, blast wave diffraction by a ramp and the one-dimensional 
shock tube problem solved on a two-dimensional grid. 

The patched-grid scheme as developed in reference 6 can be used with first- 
order accurate, explicit integration schemes. However, first-order accurate inte- 
gration schemes are insufficient to produce accurate results for a general class of 
problems. In reference 7 the patched-grid scheme of reference 6 is extended so that 
it can be used in conjunction with second-order- accurate, implicit, integration 
schemes such as the Beam-Warming scheme (ref. 5) and an implicit form of the Osher 
scheme (ref. 8). Both these integration schemes use approximate factorization to 
retain the block tridiagonal nature of the implicit equations in two and three 
spatial dimensions. The implicit patch-boundary condition, however, is developed so 
that it can be used with both factored and unfactored schemes. In reference 9, the 
explicit patched-grid scheme is extended to work with the implicit relaxation 
schemes of reference 10. Reference 11 presents two new conservative, patch-boundary 
schemes for finite-volume methods and demonstrates linear stability of these 
schemes. In reference 12 the schemes developed in reference 11 are applied in two 
and three dimensions to provide local grid refinement for airfoil and wing 
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computations. Reference 13 presents a finite-volume, patched-grid scheme for the 
Euler equations and some local grid refinement studies. Reference 14 also presents 
grid refinement studies with patched and overlaid meshes using the CSCM method 
(ref. 15). 

Reference 16 presents a conservative zonal boundary condition for problems in 
one spatial dimension. The primary emphasis of this study is the refinement of the 
mesh in regions of large solution error. In addition to refinement of the grid in 
space, the problem of using different time steps in different regions is 
addressed. A general procedure to derive consistent, conservative, patched, and 
overlaid boundary conditions is developed in reference 17. 

While the multiple-grid publications referenced above present schemes and 
results for the Euler and Navier-stokes equations, related work has also been done 
using the transonic full potential equation. Typical of this effort is the work of 
Atta (ref. 18) and Atta and Vadyak (ref. 19). The overlaid-grid approach is used in 
references 18 and 19. Reference 20 gives results obtained on overlaid grids in 
conjunction with the stream function approach. References 21 and 22 give patched 
grid results for the potential and full potential equations; the problem of conser- 
vation of fluxes is not addressed. 

In this set of lecture notes, the patch-boundary scheme of references 6, 7, 
and 9 are described in detail. The integration methods used to update the interior 
grid points are also discussed. A brief mention is made of the work on stability of 
patch-boundary schemes and the use of these schemes in Navier-Stokes calculations. 
Results are presented for inviscid, supersonic flow over a cylinder, blast wave 
diffraction by ramp, and the motion of a vortex in a freestream. These test cases 
demonstrate the quality of solutions possible with the patch-boundary scheme under 
study. 

One of the most important applications of the patched-grid approach is the 
treatment of flow regions associated with bodies which have some parts moving rela- 
tive to others, e.g., the helicopter rotor-fuselage combination, or the rotor-stator 
configurations found in turbines and compressors. It is impractical to have a 
single grid that envelops both the moving and stationary parts of the system. This 
problem can be overcome simply by containing the stationary parts in one patch (or 
set of patches) and the moving parts in another patch (or set of patches). The 
patches that contain the moving parts can be made -to be stationary relative to the 
moving parts. This approach gives rise to patch boundaries (where the moving and 
stationary patches meet) at which one set of grid points move relative to others. 

The patch-boundary treatment in this case has to be time accurate (in addition to 
the usual requirement of spatial accuracy) in order to yield reasonable results. 

The time accuracy that is possible in such patched-grid calculations is demonstrated 
iri reference 7. The test calculation of this reference has been included in this 
set of notes to demonstrate the feasibility of performing calculations on patches 
that move relative to each other. An application of this technology is also pre- 
sented. The application consists of two parabolic-arc airfoils that move relative 
to each other much like the blades in a rotor-stator combination. Results in the 
form of surface pressure histories and time varying pressure contours are included. 



As a final demonstration of the applicability of the moving-patch technique to 
realistic rotor-stator problems, it has been used to simulate the flow in a rotor- 
stator configuragion of an axial turbine (ref. 23). The numerically obtained 
results are compared with experimental data. A good agreement is obtained for both 
time-averaged surface pressure data and surface pressure amplitudes. 


INTEGRATION SCHEMES 


The patch-boundary scheme is developed for explicit first-order accurate inte- 
gration schemes and for implicit, relaxation schemes that are both first- and 
second-order accurate in a later section. In this section, the explicit and 
implicit versions of the Osher scheme (refs. 8, 10, and 24) are presented. Although 
the flux-linearizations that are required for the implicit technique are developed 
with the Osher scheme in mind, they can be used equally effectively with the split- 
flux scheme and Roe's scheme (ref. 25). The implicit scheme that is presented is 
iterative in nature, but reverts to a conventional noniterative implicit scheme when 
the number of iterations is restricted to one. The iterative procedure considerably 
simplifies the extension of explicit patch-boundary schemes to implicit patch- 
boundary schemes (details of which are given in the next section). 

Consider the unsteady Euler equations in two dimensions 


Q, + E + F =0 
t x y 


(1) 


where the vectors Q, E, and F are given by 

fp \ / pu 


Q = 



E = 


P + pu 
puv 


F = 



(2) 


^(e + p)uy 

where p is the density, p is the pressure, u and v are the velocities in the x 
and y directions respectively, and e is the total energy per unit volume 


e = 



(3) 


Establishing the independent variable transformation 

t = t 

£ = £(x,y,t) 
n = n(x,y,t) 


(4) 
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and applying this transformation to equation (1), the following is obtained 

0 + E„ + F =0 
t 5 n 

where 


( 5 ) 


Q = Q/J ^ 

E(Q,C) = (C t Q + 5 X E + £ F)/J 


F(Q,n) = (n,Q + n E + n F)/J 

u a y 

J : 5 1 - I] 5 

s x y x^y 


J 


( 6 ) 


The notation E(Q,^) and F(Q,n) is used to show the dependence of these quantities 
on the metrics of the transformation. 

A conservative finite-difference scheme for equation (5) is given by 


- F m 


AH n P m R m R m 

liiV2^k ~ ^,1-1/2 ,k r ,1,k+1/2 ~ M.k-1/2 
Ax + H An = 


(7) 


where E. . . and F. . . /0 are numerical fluxes consistent with the transformed 

J+1/2,k j,k+1/2 

physical fluxes £ and F. The difference scheme (eq. (7)) is explicit when m = n 
and fully implicit when m = n + 1 . 


The First-Order Accurate Osher Scheme 
The numerical flux for the first-order Osher scheme is given by (ref. 24) 


1 


E j+1/2,k = 2 [E(Q j,k’ 5 j + 1/2,k ) + E(Q j + 1,k’ 5 j + 1/2,k 


where 


AE 


- AE + (Q jjk ,Q j+ 1 >k ,C J+1/2>k ) + AE (Q J>k ,Q j+ 1 )k »C j+1/2)k 


± r Q j-Hl ,kr 3 g -It 

(Q j,k ,Q j+1 ,k ,5 j+1/2,k ) = J [aQ (Q,5 j+1/2,k ) J dQ 


)] ( 8 ) 


(9) 


Q j » k 
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The numerical flux F. . . can be obtained in a similar manner. The evaluation 

j ,k+1/2 

of the line integral can be found in reference 24. 

An evaluation of the numerical flux in equation (8) using the dependent 
variables at the nth level results in an explicit scheme and difference equations 
that are linear. However, an evaluation of these fluxes at the (n+1)th time level 
results in an implicit scheme and difference equations that are nonlinear and need 
to be solved in an iterative manner. The usual strategy that is employed at this 
stage is the linearization of the numerical fluxes with respect to the time like 
variable x. The resulting system of linear equations is then solved in order to 
update the dependent variables. The linearization process depends on the scheme 
used to determine the numerical fluxes. The linearization of the Osher fluxes is 
discussed below. 


The time linearization of the numerical flux E n+ ] . requires the lineariza- 

j+1/2,k 

tion of the integral given in equation (9). The linearization of this integral is a 
cumbersome process that is computationally expensive. Hence, the following approxi- 
mation (ref. 8) is made to simplify the linearization process: 


[aQ (Q,5 j+1/2,k ) ] = 8Q (Q ’^j+1/2,k )] 

where 

E ( Q ,5j + i /2 , k ) = [jQ (Q ’ 5 j+1/2,k ) J Q 


(10) 


(ID 


Substituting equation (10) into equation (9) yields 


UE X*’W'h*i/^ in+1 = <e k* - 


( 12 ) 


Equations (8) and (12) together yield 


(E j+1/2,k ) approx = ^ (Q 


j , k’ ^ j+1/2 , k } + E (Q j + 1,k’ C j+1/2,k )J 


n+ 


(13) 


Linearizing equation (13) with respect to x and making use of equation (10) once 
again, it can be shown that 


<E n *' 


j+1/2,k approx 


= (E n ) 

j+1/2,k approx 


< 5 12 


n 


* 


(A 


j+1 »k 


) n aQ, 


J+1 .k 


(14) 
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where 


5 J’ k = Lo (Q j,k’ 5 j+1/2,k ) ] (15 

The numerical flux for the first-order accurate Osher scheme using this approximate 
linearization can now be written as 


( R n+ ^ ^ 

j+1/2,k'linearized ” j+1/2,k 


[S (Q j,k’ 5 j+1/2,k )A ^j,k 


S (Q j+1 ,k’^j+1/2,k )A ^j+1 ,k ] 


( 16 ) 


where E n is evaluated using equation (8). It should be noted that the use 

J+ 1 / 2 , k 

of the approximate linearization does not result in any loss of the conservative 
properties of the scheme. 


The iterative implicit technique of reference 8 as applied to the first-order 
accurate Osher scheme now takes the form 


I + ~ (v r A + . + A r. . ) 
A? K J,k K J,k 


y At 


= -Ax 


j 1 (v §* . 

An n j,k 

?P 


1 n B j 




Q n E 

+ b >1/2, k 




rP rP rP \ 

+ j,k + l/2 ' b J,k-1/2 ^ 


An 


(17) 


where Q p is an approximation to Q n+1 . When p = 0, Q p = Q n and when equa- 
tion (17) is iterated to convergence at a given time-step, Q P = 0 n+ \ It should be 
noted that, because the left-hand side of this equation can be made equal to zero at 
each time step (by iterating to convergence), linearization errors can be driven to 
zero during the iterative process. For problems where only the asymptotic steady- 
state is of interest, the iteration process need not be carried to convergence at 
each time-step. In fact, when the number of iterations is restricted to one, the 
scheme reverts to a conventional noniterative scheme of the type in reference 5 (but 
unfactored) . 


Unfortunately, equation (17) is extremely time consuming to solve in a direct 
fashion because of the large bandwidth of the matrix on the left-hand side. At this 
point, two options are available to the user; the first option is to use approximate 
factorization as in references 5 and 8, the second option is to use a relaxation 
strategy as in references 9 and 10. The approximately factored form of equa- 
tion (17) is given below 
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( 18 ) 


At 

A 5 


(v r A . 

5 j,k 


A A" )"1 P fl + (V B* . 

5 j,kj L An n J,i 



= RHS of equation (17) 

Clearly, approximate factorization reduces the bandwidth (the single large band- 
width matrix is converted into two block-tridiagonal matrices). If the iteration 
process is carried to convergence, the factorization error is driven to zero. 
However, at very large time steps, factorization error may cause the iteration 
process to diverge or enter a limit cycle. 

Relaxation schemes on the other hand do not require factorization and hence, 
typically, permit the use of much larger time steps. The iterative form of equa- 
tion (17) is ideally suited to relaxation schemes. Although there are a host of 
different relaxation schemes available, the ones that have been chosen for this 
study are schemes that are most amenable to vectorization ; in particular, pointwise 
and line relaxation schemes of the non-Gauss-Seidel type. The pointwise relaxation 
scheme is obtained by discarding all the nondiagonal terms on the left-hand side of 
equation (17). This yields 


I + 


At 

ae; 




- a 7 . ) 

J,k 


At /s+ 

An j ,k 


r 


(Q 


P+1 aP 
J,k " 


) = RHS of equation (17) 


(19) 


The line relaxation scheme is obtained by retaining all the diagonal terms and terms 
corresponding to points on a given line, for example, a constant n line. This 
yields 


I + 


£ ( v St 

AE. K j,k 


A A. . ) + 

5 j,k' 


At 

An 






= RHS of equation (17) (20) 

If the constant ? lines had been chosen instead, the resulting equations would be 


fl + ~ (A* - At . ) + (V B't 

L ae. j,k j ,k An n j,l 


A n B j,k>J 


nP +1 


(Qt 


j »k 



) 


= RHS of equation (17) (21) 

Variants of the line and pointwise relaxation schemes (such as the zebra and check- 
erboard schemes) are also discussed in reference 10. 


9 



Although equation (17) is fully conservative (in spite of the approximate 
linearization), equations ( 19)— (21 ) are not conservative in time unless they are 
iterated to convergence at each time step (this would be necessary for time depen- 
dent problems where the transients are of interest). However, for problems where 
only the time asymptotic solutions are required, equations (18)~(20) yield solutions 
that satisfy the conservation laws at steady-state without requiring iteration to 
convergence at each time step. This is because the right-hand side of these equa- 
tions is fully conservative. 


The Second-Order-Accurate Osher Scheme 

The numerical flux for the second-order accurate Osher scheme is given by 
(ref. 24) 

E, , , = E . , , (first order Osher scheme) 

j+1/2,k j+1/2,k 


[AE + (Q._ 


1,k jQ j,k’ C j + 1/2,k ) “ 


AE (Q j + 1,k’ Q j4-2,k’ 5 j+1/2,k )] (22) 


where AE + are evaluated as before. Mote that the numerical flux for the second- 
order scheme is denoted by E^ + ^ 2 ^ to distinguish it from the numerical flux for 

the first-order Osher scheme. Linearization of all the terms in equation (22) would 
result in block-pentadiagonal matrices for line-relaxation schemes. Hence, only the 
terms corresponding to the first-order scheme and those second-order terms that 
contribute to the diagonal elements of the matrix are linearized (private communica- 
tion, S. R. Chakravarthy , Rockwell International Science Center). The resulting 
iterative, implicit scheme takes the form 
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(i + Al ri (A 4 : . - A" ) + V a + . + A r . 1 + Al [1 (g + _ g- ) + 7 B + 

l as L2 j,k j,k s j,k 5 j,kJ An l_2 j,k J,k n j,k 

/qP _ Q n gP _ gP 

nP x At [ i’ k w 1,k b j+1/2,k j- 1 /2 , k 

“ y j,k ; - “ Ax \ At AS 

k-1/2 ^ 

- — [AE (Q 3 . 1 >k »Qj fk »Cj +1/2|k ) - AE (Q j- 2 ,k ,Q j- 1 ,k’^j- 1 / 2 ,k )] 


f I . ft f 


\ AS L 

+ A B. . j 

} P <«' 

n J,KJ 

/v 

pP 


M.k+1/2 

T 

An 


At 


2AS [AE ( Q j+1 , k’ Q j+2 , k’ ^ j+1 /2 , k } ' AE (Q j , k’ Q j+1 , k’ ^ j-1 /2 , k 


)] 


' 2An UF (Q j,k-r Q j,k’ n j,k+1/2 ) “ AF (Q j ,k~2’ Q j , k-1’ n j , k- 1 /2 } 


At - - 1 

+ 2M [AF (Q j,k+1 ,Q j,k+2’ r| j,k+1/2 ) " AF (Q j , k’ Q j , k+1 ’ n j , k-1/2^ 


(23) 


The terms that result in second-order spatial accuracy (on the right-hand side of 
eq. (23)) have been modified, that is, these flux differences are obtained from a 
flux-limiting process. Details of flux-limiting for the Osher scheme can be found 
in reference 24. Flux-limiting is essential to provide the diagonal dominance 
required for relaxation schemes (ref. 10). 

The pointwise relaxation scheme is now obtained by discarding all nondiagonal 
terms, that is, 


I 


+ 


Mi a + 

2AS 1 J,k 








= RHS of equation (23) 

and the line relaxation scheme corresponding to equation (20) takes the form 

P 


(24) 


XI + 


45 [2 <A j,k: ~ fi j,k ) *■ Vj,k * 4 5 5 j,k ) ] + 24 n <B j,k ‘ 8 J , k ^ 




= RHS of equation (23) 


(25) 


The comments regarding the conservative properties of the first-order schemes apply 
to the second-order schemes as well. As stated earlier, although equa- 
tions ( 1 7 ) — ( 25 ) have been developed for the Osher scheme, they are equally 
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applicable to the split-flux scheme and Roe's scheme; only the evaluation of the 
fluxes on the right hand side of these equations changes according to the particular 
integration scheme chosen. 


THE PATCH BOUNDARY SCHEME 


In order that information is transferred accurately and easily from patch to 
patch, a patch-boundary condition must satisfy several requirements. Some of these 
are listed below.. The patch-boundary condition must be: 

1 . Numerically stable 

2. Spatially and temporally accurate 

3. Conservative so that flow discontinuities can move from one grid to another 
without any distortion 

4. Easily applicable in generalized coordinates. 

A patch-boundary condition that satisfies these requirements and that can be 
used with explicit and implicit schemes is presented in this section. 


The Explicit Patch-Boundary Scheme 

In order to develop the explicit patch-boundary condition, we consider the 
unsteady Euler equations in two spatial dimensions (eq. (1)). Equation (7) presents 
a conservative-difference scheme for the unsteady Euler equations. For the simple 
transformation of coordinates 5 = x, n = y, and x = t, this set of equations 
(eq. (7)) can be written as 

AQ j,k t E j+1/2,k ~ E j-1/2,k , F j , k+1 /2 ~ F j,k-1/2 . Q 6 

At Ax Ay v ' 

/V A 

where E and F are evaluated at the nth time step. Equation (26) can be alterna- 
tively written as 


( Q n+ 1 _ Q n ) + 

At ^j,k y j,k ; + 


Ay(E 


J+1/2,k 


E j-1/2,k } 


+ 4X(F j,k.t/2 


- F ) = 

j,k-1/2 ; 


(27) 


In this form the various terms in the equation can be easily interpreted; for exam- 
ple, the first term 
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in the cell ABCD (fig. 5). 


AxAy 

At 


<ft 



represents the rate of increase of the variable Qj k 
The term 


-Ay(E 


j+1/2,k 


E ) 

j-1/2,k' 


represents the net influx of mass, momentum, and energy into the cell through the 
sides AB and CD, and the term 

A /V 

~ Ax(F j,k+1/2 " F j,k-1/2 ) 


represents the net influx of the same quantities through the sides AC and BD. 
A summation of the term 


AxAy 

At 


‘ft - ft> 


over all the grid points yields (using eq. (27)) 


jmax kmax 


_ y y ami (Q n+i _ n n ) 

Lj At ^ y j,k w j,k ; 


J=1 k=1 

jmax 

4x Z Z, 

j=i 


kmax 

1/2 ' F j,kmax+1/2 ) + Ay 2 (1 

k=1 


1/2, k _ E jmax+1/2,k ) 


( 28 ) 


that is, S is only a sum of the boundary fluxes; the interior fluxes cancel each 
other out because the difference scheme (eq. (26)) is conservative. Equation (28) 
represents the global conservation property of any scheme that can be reprerented as 
in equation (26). 


Consider the grid shown in figure 6. The line AB represents the patch boundary 
that separates the two grids that are used to discretize the given region. Let 2, 
and m be the indices used in the x and y directions, respectively, in patch 1 and 
let j and k be the corresponding indices for patch 2. Let n represent the time 
step for both patches. A superscript within parentheses will denote the patch to 
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which a given quantity belongs, e.g., Ax^ ^ denotes the mesh spacing in the x 
direction in patch 1. 

One condition that must be satisfied across the patch boundary is the continu- 
ity of the dependent variables. This condition can be easily satisfied by integrat- 
ing the equations of motion to update Q on one side of the patch boundary and 
interpolating these variables to obtain the updated variables on the other side of 
the patch boundary. This interpolation results in the dependent variables being 
continuous (across the patch boundary) to the order of accuracy of the interpolation 
scheme. A linear interpolation scheme was used to obtain the results presented in 
this set of notes. 


Assume that the patch-boundary points of patch 2 are to be updated using the 
finite-difference scheme of equation (26). This calculation requires the fluxes 
"( 2 ) 

Fj The Question of how to calculate these fluxes introduces the second condi- 

tion that is to be satisfied at the patch boundary in a natural way; that is, these 


fluxes have to be calculated such that global conservation is maintained. In other 
words, a summation of (Ax^ Ay^/At)AQ^ ^ and (Ax^^Ay^/At)AQ^ over all the 


&,m 


J.k 


cells in the region of interest (as was done in eq. (28)) should once again result 
in only the boundary fluxes. The sum should not contain any residual fluxes near 
the patch boundary. Mote that cells corresponding to the patch-boundary points of 
patch 1 are not to be included in this summation. This is because the area that 
these cells represent has already been accounted for with the inclusion of the cells 
corresponding to the patch-boundary points of patch 2 (a sum of the cell areas 
included in the above summation should result in the total area covered by the two 
patches). A typical cell (RSTU) of a patch boundary point (j,1) is shown in 
figure 6. The points R and S are midpoints of the cells they lie in while the 
points T and U are obtained as follows: The constant j lines of patch 2 are 
extrapolated into patch 1 to intersect the line CD (CD corresponds to 
m = mmax-1/2 in patch 1 and to k = 1/2 in patch 2). The intersection points have 
the indices (J, 1/2) . Point T is midway between the points (j + 1,1/2) and (j,1/2) 
while point U is midway between points (j,1/2) and (j - 1,1/2). 


The global conservation property can be shown to be satisfied if the following 
relationship is satisfied 


. ( 2 ) 

Ax / 17 ( 2 ) 


;( 2 ) 


(f; ' ._ + F . , .1 + Ax 

1,1/2 jmax,1/2 


jmax-1 

(2) z 

j=2 


( p(i) 


+ f! 1} 


n(2) 

F j,1/2 


2max-1 


2 1, mmax-1/2 8,max, mmax-1/2 


,) + Ax 


(D 


s 

a=2 


;(D 

2., mmax-1/2 


(29) 
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A close examination of equation (29) shows that each side of this equation is 


nothing but a discrete form of the line integral of the numerical flux F along the 
line CD in figure 6 while the equation itself represents flux conservation across 
the patch boundary. Equation (29) is only a necessary condition and is not 


sufficient to define the fluxes 


F< M/2 canno &,£f 

specify the F. ' 

J i ' ' c- 


A (2) 

Fj '^2 in a physically meaningful fashion (the 


obtained from eq. (29) alone since eq. (29) does not uniquely 

). 


Assume that the 


"( 2 ) "( 1 ) 

Fj j /2 are obtained by interpolating the Pj >I J imax _ 1/ 2 » i-e*» 


p(2) 

j, 1/2 


■ E 


"( 1 ) 

N F v 

j , 2. «,,mmax-1/2 


(30) 


sup 


where the N, ^ are interpolation coefficients and p and q define the set of 
fluxes of patch 1 that will be used in the interpolation. We now describe a very 
simple way of obtaining the interpolation coefficients Wj ^ such that equa- 
tion (29) is automatically satisfied. Let the line CD in figure 7 correspond to the 
line CD in figure 6. The dots represent the grid points of patch 1 and the crosses 

A ( 1 ) 

represent those of patch 2. Representative numerical values of F' . are 

v £,mmax-1/2 


plotted on the positive y axis. 


Assume a piece-wise constant variation of 

.( 1 ) 


A ( 1 ) A ( i ) 

F , that is, F; , is constant between x' 1 ( and 

2,,mmax-1/2’ ’ !l,mmax-1/2 a-1/2,mmax-1/2 

( 1 ) 

x ?,+ 1/2 mm x 1/2' Consider a patch-boundary point of patch 2, (j,1). Let E 
midway 'Between (j - 1,1/2) and (j,1/2) and F be midway between (j,1/2) and 

"( 2 ) 

(j + 1,1/2). The Fj are now calculated from 


be 


Ax 


( 2 )"( 2 ) 
j, 1/2 



p(1) 

a,mmax-1/2 


dx 


or 


:(2) 

j, 1/2 


Ax 


(2) 


r 


pd) 

8,,mmax- 


1/2 


dx 




pd) 

J,,mmax-1/2 


«,=p 


(3D 
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where the values of Nj ^ are given by 


M. = { 

J.a 


0 if x (1) x (1) < x (2) 

4+1/2' ft- 1/2 ~ X j — 1/2 

0 if x (1) x (1) > x (2) 

X ft+1/2’ X ft-1/2 ~ j+1/2 

(x r " X ft )/(x >1/2 " X j?i/2 ) otherwise 


( 32 ) 


and 


x - min(x^ 2 ^ x^ ^ ) 

r " min X j+1/2 ’ x ft+1 /2 


h s max(x j-!/2' x ^i/2 > 


The simple expressions for N, 0 given in equation (32) are valid only for a piece- 
wise constant variation of the numerical fluxes. A piecewise linear or any other 


variation would result in different formulae for the 


j»*: 


Equation (32), in an 


indirect manner, also yields the end points in the interpolation, p and q (p and q 
only include that set of fluxes of patch 1 that are multiplied by nonzero interpola- 

;( 2 ) 


tion coefficients for a given flux of patch 2 ), 
are calculated as 


The end fluxes 


p<2> 

jmax, 1/2 


;(2) 

1 , 1/2 


Ax 


(2) 


F 1 , 1/2 a " d 


.0 


F (,) dx 

ft,mmax-1/2 


and 


E 


ft =1 


M ; 

1 , ft ft,mmax-1/2 


(33) 


;(2) 

jmax, 1/2 



"(1 ) 

F dx 

ft,mmax-1/2 ax 


ftmax 


■ E 


N 

jmax, ft ft,mmax- 


1/2 


(34) 


ft=u 
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The fluxes thus obtained satisfy equation (29) automatically. The shaded areas in 
figure 7 represent the values of the integrals in equations (31), (33), and (34). 

Summarizing, the explicit patched--grid scheme consists of the following three 
steps : 

1. Integrate the dependent variables at grid points (of both of the grids) 
that do not belong to the patch boundary. 

2. Integrate the dependent variables at the patch-boundary points of one of 
the patches (for instance, patch 2 of fig. 6) using a scheme that conserves fluxes 
across the patch boundary. 

3. Obtain the dependent variables at the patch-boundary points of the other 
patch (for instance, patch 1 of fig. 6) such that the dependent variables are con- 
tinuous along the patch boundary. 

The foregoing discussion pertained to simple rectangular grids. The extension 
of the method to arbitrary curvilinear grid systems is straightforward and is out- 
lined below. Consider the curvilinear grids used to discretize the region shown in 
figure 8. Establishing two independent variable transformations 

= t 

e; ( i) = s (l) (x,y,t) (35) 

= r/ 1 ^ ( x , y , t ) 

and applying these transformations to equation (1) we obtain 


+ E ( J-\ + = 0 

x (l) 5 (1) n 1 


(36) 


where 


0 (l) = Q/J (l) 


E (l) (Q,S (l) ) 


F (l) (Q,n (l) ) 


+ 4 i} e + 4 i W i} 


(n t ( . l)Q + 4 1} E + n^ l) F)/J (l) 


(37) 


,( i) _(i) (i) ( i)_(i) 

’ - ^x n y " n x 5 y 
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The notation E ^ ^ ) and F * (Q,n ^ ) is used to show the dependence of these 

quantities on the metrics of the transformation. Let the conservative difference 
schemes used to integrate equation (36) be given by 


AQ 


(1) ;d) 


Sljn 


( 1 ) 


At 


( 1 ) 


J g.+ 1/2,m 8, -1 /2 ,m 


A 5 


( 1 ) 


:d) ;o) 

«,,m+1/2 a,m~1/2 


An 


(1) 


(38) 


and 


a0 ( . 2) 


g(2) 


lUxk + j + 1/2,k 

/ o \ ^ 


At 


(2) 


A5 


( 2 ) 


p(2) ?(2) p(2) 

JzJ/fLk A ^kH-1/2." j,k-1/2 = 0 


An 


(2) 


(39) 


The interior points of each patch are updated using the appropriate metrics and 
dependent variables. At the patch boundary, once again, the grid points of one 
patch are updated by integrating the equations of motion (for instance patch 2) and 
those in the other patch (for instance patch 1) are updated by interpolating the 
dependent variables of the first patch. 


A typical cell RSTU associated with the patch-boundary point (j,1) of patch 2 
is shown in figure 8. The points R, S, T, and U are defined as in the previous 
case. The metrics of the transformation at the point (j,1) are defined in a manner 
consistent with the shape and size of the cell RSTU, that is, 


U (2) > j , 1 
n 


(y „(2) ) j,1 


<X e <2) ) j , 1 




x (2) - x <2) 

j ,3/2 X j,1/2 


(2) (2) 
y J,3/2 " y J,1/2 

1 / (2) (2) . 
2 ^ X j+1 , 1 " X j~1 , 1 ^ 


1 fv (2) v (2) 

2 {y j+1,1 " 


(40) 


) 


( 2 ) 1 . ( 2 ) ( 2 ), 

X j,3/2 = 2 X j ,2 + X j,i } 


/2) 
j ,3/2 


•1 fv (2) , v (2) ) 

2 (y j,2 + y j,1 } 


• Flux conservation across the line CD (as in the previous case) requires that 
the following condition be satisfied: 
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jmax-1 


1 /p(2) . p(2) \ + 

2 (F 1,1/2 jmax, 1/2 


y f (2) 
j > i/2 


j=2 


Slmax-1 


= 1 (F (1) + p(D 

2 1,mmax-1/2 Slmax, mmax- 1/2 


) + 


Z f *’ 


) 

mmax- 1 /2 


(41) 


a -2 


Equation (41) assumes that A^ 1 = 1 and An 1 =1. In order to satisfy equa- 
tion (41), a running parameter s is established along the line CD. The quantity 
s respresents the distance of a point from the point C along the curve CD. 

Figure 7 shows the line CD stretched into a straight line along the s axis. 

;d) 

!,mmax-1/2 

^ is assumed between grid 

1, mmax- 1/2 & 


Representative numerical values of F) 
axis and a piecewise constant variation of F) 


are plotted on the positive y 


points as before. Defining the points E and F as before, equation (41) can be 

"( 2 ) 

satisfied exactly by evaluating F^ ' /2 from 


p ( 2 ) , 0 ( 2 ) , . 

F j , 1/2 Q,n j , 1/2 “ 


f 


p(D (0 n (1) Ids 

1, mmax- 1/2^ ’ n 2,,mmax-1/2 ; 

s (1) - s (1) 

1 + 1/2 1 - 1/2 


I 


a=p 


W F (1) (0 n (1) 1 

j,l l,mmax-1/2 v ’ 1, mmax- 1/2' 


where the values of Nj ^ are given by 


(42) 




■o if s (,) ,<') < s (2) 

1 + 1 / 2 ’ 1 - 1/2 “ j - 1/2 

0 if s (,) s (,) > s <2) 

U 1 11+1/2’ «-1/2 " j + 1/2 

(s r - V /( 41l/2 - 4)!/2> otherwise 


(43) 
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and 


s - min(s^ ^ ) 

r ' { j+1/2’ i+]/2 > 


2 ) ( 1 ) \ 
- max(s J _ 1/2 ,s Jl _ 1/2 ) 


( 1 ) ( 1 ) 

The term (s ' /0 - 0 ) in the denominator of the integral in equation (42) 

I / <-. I / a / ^ \ 

serves to convert the numerical flux F„ . into a flux per unit length (the 

2.,mmax-1/2 ^ 0 

metrics contained within the numerical flux take the length of the side of the cell 
into account). The integration process reintroduces the length parameter. 

-'( 2 ) 

The preceding method of obtaining the F. . (eq. (42)) is not freestream pre- 

J > ' ' £ - 

serving. However, a repeated integration of the governing equations with freestream 
conditions everywhere (as initial conditions), and without boundary conditions on a 
grid of the type shown in figure 8, resulted in only a 0.1% drift in the density in 
the vicinity of the patch boundary. It can be shown (ref. 6) that the drift in 
freestream conditions near the patch boundary is proportional to the curvature of 
the patch boundary, and is caused by terms that are second order in magnitude. 

Hence, it behooves the user to use patch boundaries with moderate curvature. 


The patch-boundary technique has thus far been developed for a region divided 
into two adjacent patches. The technique, however, is not restricted to two-patch 
calculations, but can be generalized to an arbitrary number of patches positioned 
arbitrarily with respect to each other. As an example of using multiple patches, 

.the region shown in figure 9 is divided into four patches with the use of three 
patch boundaries. The lines along which a flux balance is carried out are AC, FG, 
DC, and EB. The shaded areas depict typical cell areas for points situated on the 
patch boundaries. Grid points that belong to only one patch-boundary usually have 
four-sided cell areas. However, patch-boundary grid points that are close to the 
intersection of two or more patch boundaries may sometimes be required to have cell 
areas that are bounded by more than four sides (point I in fig. 9). These nonstan- 
dard cell shapes are required to cover the entire region by cell areas corresponding 
to grid points which are updated by integrating the equations of motion. The cell 
areas corresponding to patch-boundary points can be easily determined once the lines 
along which a flux balance is carried out are established. The evaluation of the 
fluxes necessary to integrate the nonstandard cells is very similar to the procedure 
used for standard four-sided cells. Details regarding the treatment of nonstandard 
cells can be found in reference 6. 


"( 1 ) 

The foregoing discussion assumed that the numerical fluxes F; . were 

& 6 £,mmax-1/2 

readily available. For the first-order scheme an examination of equation (8) will 

( 1 ) ( 1 ) ( 1 ) 

show that only the quantities Q; . , Q; , and n„ , are sufficient 

M il,mmax-1’ 8,,mmax’ 8,,mmax-1/2 

to define F„ . /0 . Since these quantities are easily available, calculation 

S,,mmax-1/2 J ’ 


20 



of F 


( 1 ) 


is a simple matter. However, for the second-order accurate scheme, 


5,,mmax-1/2 

equation (22) shows that the flux difference 


AF"(Q 


(1) 

A,mmax 


,Q 


(D 

2,, mmax+1 


(1) 

’ n 2,,mmax-1/2 


) 


is required to calculate F; . The calculation of this flux difference, in 

M S,,mmax-1/2 ’ 

( 1 ) 

turn, requires the vector Q; . (the line mmax+1 in patch 1 corresponds to 

’ M a, mmax+1 (d 

the line k = 2" in patch 2). The vector Cr . is determined simply by 

% ^ mmax+ i 

extrapolating the constant % lines into patch 2 to intersect the line k = 2 in 

,(D 


patch 2. The values of Q 
\( 2 ) 


mmax+1 


are then interpolated from the values of 


Q. ' . A simple linear interpolation was used to obtain the results presented later 

J i *- 

in these notes. An alternative to this procedure, and one that avoids the interpo- 
lation procedure described above, is to calculate the flux differences 
( 2 ) ( 2 ) ( 2 ) 

AF (Q. . ,Q . „,n. 1/P ) and to perform a flux balance as follows: 

J y ‘ J t J 9 I 


AF"(Q 


( 1 ) 


8,,mmax 


,Q 


(D 

St ,mmax+ 


( 1 ) X 

* n S,,mmax-1/2 ; 






1 /2 


(44) 


J=P 


Notice that in equation (44) the flux balance is being performed in the opposite 
direction (compared to the flux balance described earlier). Similar procedures, 
that is, interpolation or an additional flux balance, are required to calculate the 
"( 2 ) 

numerical flux F. also. Both the methods, which are described above, do not 

J } d/<- 

in any way affect the conservative property of the patch-boundary scheme. The 
additional flux difference terms that are required to make the scheme second-order 
accurate are once again perfectly balanced across the patch boundary. 


The Implicit, Relaxation Patch-Boundary Scheme 

Now that the calculation of the boundary fluxes has been described, all that 
remains to be seen is how these ideas fit into the implicit, relaxation schemes 
given by equations ( 19) — (21 ) and ( 24 ) - ( 25 ) - The right-hand side of all these equa- 
tions can be calculated at all grid points in the interior and at the patch bound- 
ary, using either the usual definition of numerical fluxes or in the special manner 
described in the previous subsection. 

The pointwise relaxation schemes given by equations (19) and (24) remain com- 
pletely unchanged except for the right-hand side at patch boundaries. This is 
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because spatial derivatives of the matrices A and B do not appear on the left- 
hand side of these equations. For time-asymptotic problems, the number of itera- 
tions per time step can be restricted to just one. If for reasons of computational 
efficiency (increased convergence rates) it becomes necessary to use more than one 
iteration per time step, then the most reasonable approach seems to be to perform 
each iteration in all the patches before performing the next iteration. 

For the example shown in figure 8, where the constant 5 lines are discontin- 
uous at the patch boundary, the line relaxation schemes given by equations (20) 
and 25 also remain unchanged as in the case of the pointwise relaxation schemes. 

This is becausethe n derivatives of the matrices B + and B~ do not appear on the 
left-hand side of these equations. However, the line relaxation scheme given by 
equation (21) needs to be modified at the patch boundary of patch 2 because the 
backward difference V B. . is not defined. Several options are available to 

n J,k F 

overcome this difficulty. The simplest solution would be to revert to the pointwise 
relaxation scheme at the patch-boundary points, that is, 


i * it - *;,.<> - tt 


(qP + ^ -#?.)= RHS of equation (17) (45) 

J 1 K J J K 


when k = 1 , and 


1 * ft - 5 L*> * ft "tf,* 


W]' 


8 j,k> 


= RHS of equation (17) (46) 


when k / 1 . A second approach would be to retain only the diagonal terms from the 

derivative V St . which then yields 

n j,k 

[ * ft (J ),k - S L*> * ft (8 f,k * Vfypjft - 8 j,*> * RHS of equatlon <17) 


(47) 


when k = 1 , and 


T 1 + 77 - A" ) + I 1 (7 S + . + A ST „)l P (Q l Tl - 3 P , 

L J> k n j,k n j,kj j,k j,l 


= RHS of equation (17) (48) 


when k i 1. The results presented in this study (and obtained with the line relax- 
ation schemes) were calculated using equations (47) and (48). The comments made 
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earlier regarding the number of iterations per time step (for the pointwise 
relaxation schemes) are valid for the line relaxation schemes also. 

The implicit, relaxation, patch-boundary scheme can be summarized in the fol- 
lowing steps: 

1 . Integrate the dependent variables at all the grid points of patch 2 using 
the pointwise or line relaxation schemes with their special forms at the patch 
boundary (only one iteration). 

2. Interpolate the newly obtained values of [(q! 2 :!) P+ ^ - (Q^i!)^] along the 

patch boundary to yield a new set of values of [ (q'P.L,.) 134 ^ - (q!^__„)^]. 

56 j mma x 56 j mm a x 

3. Integrate the dependent variables at the grid points of patch 1 using the 
pointwise or line relaxation schemes (only one iteration) and the most recent values 

°f [(Q^mmax ^ 1 " ^a]mmax^ Dirichlefc “ t yP e boundary condition). 

4. Interpolate the values of (Q^. 2 ]) p+1 to obtain the values of (q[^ „ v ) P+1 

j , I 56 j mmax 

(discard the ones obtained as a result of the integration). 

5. If the maximum value of the magnitudes of all [(Q^) p+1 - (Q^) p ] is less 
than the prescribed tolerance limit, go to the next integration step; if not, go 
back to step 1 and reiterate. 

Step 2 and the implementation of the Dirichlet boundary condition in step 3 are 
not required for pointwise relaxation schemes and line relaxation schemes where the 
lines are chosen to be of the same family as the patch-boundary line. In the case 
of implicit, relaxation, patch-boundary schemes one iteration per time step has been 
found to be sufficient to maintain stability, multiple iterations being required to 
obtain time accuracy. 


The Factored, Implicit Patch-Boundary Scheme 

A linearized, implicit scheme requires the linearization of the numerical 
fluxes. The scheme then entails the solution of a system of linear equations where 

the unknowns are the incremental changes in the dependent variables at the grid 

points. The matrix associated with this system of linear equations is sparse, but 

may have a large bandwidth and is, hence, computationally expensive to solve. The 

strategies used to overcome this problem are approximate factorization and relaxa- 
tion procedures. In relaxation procedures, inter-grid point connections are altered 
using certain approximations which are then accounted for in the iteration process 
that goes hand-in-hand with relaxation procedures. The new grid point connectivi- 
ties are chosen such that the associated matrix is much simpler to solve (for exam- 
ple, the matrix is diagonal for pointwise relaxation schemes). In factored, 
implicit schemes the connectivities between various grid points is retained, but 
additional connections are brought into play (the factorization error term) so that 
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the original matrix can be expressed as the product of two or more matrices that are 
easier to solve (for example, two tridiagonal matrices in eq. (18)). 

A linearization of the numerical fluxes in the vicinity of the patch-boundaries 
results in the incremental changes in the dependent variables on the two sides of 
the patch boundary becoming coupled. The factored, implicit scheme cannot be fac- 
tored at the patch boundary. Hence, we resort to the kind of approximations made in 
the relaxation approach (only at the patch-boundary) and use the iterative technique 
to correct for the approximations made. Details regarding the factored, implicit 
scheme can be found in reference 7. 


ADDITIONAL TOPICS 


In this section, we briefly touch upon various topics that are directly related 
to patched-grid simulations. While these topics are no less important than the 
patched-boundary method itself, the scope of the current study does not permit an 
exhaustive development of each one of them. 


Stability of Patched-Grid Schemes 

The numerical stability of the patched-grid scheme developed in the previous 
section is demonstrated in a practical manner in the next section by applying it to 
several example problems. The examples include the motion of strong shocks and 
other flow discontinuities through patch boundaries (both moving and stationary 
patches) and, therefore, they amply test the patch-boundary scheme. However, a more 
formal analysis of stability (even if it is only for the patched-grid scheme applied 
to a linear differential equation) is a necessity. Reference 11 presents such an 
analysis for the patch-boundary condition presented here and also for some new ones 
developed for the finite-volume technique. 

Stability analyses for interior grid points are usually performed using a 
Fourier technique. An alternate method is required to analyze the stability of a 
scheme when the boundaries of the computational region are also to be included in 
the analysis, for example, patch boundaries. Reference 11 uses the amplification 
matrix to test the stability of the interior scheme together with the patch-boundary 
scheme. This is done by computing the largest eigenvalue and the L 2 norm of 
arbitrary powers of the amplification matrix. The scheme is stable if and only if 


| | [G( Ax) ] n | | < K e anAfc ; n = 1,2, — 

where G is the amplification matrix, Ax is a measure of the grid spacing, and K 
and a are some constants independent of Ax. Unlike traditional stability analy- 
ses, using Fourier techniques, this approach relies heavily on numerical methods to 
compute eigenvalues and L 2 norms of matrices. Reference 11 demonstrates stability 
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of the patched-grid schemes presented earlier (as applied to a system of linear 
hyperbolic partial differential equations) and also sheds light on the processes 
that govern the stability of patched-grid schemes. 


Patched-Grid Calculations in Three Dimensions 

The patched-grid scheme can be extended in a straightforward manner to three 
spatial dimensions. Although the programming logic in three dimensions is more 
complex, conceptually there are no problems in extending the two-dimensional scheme 
to more dimensions. Preliminary results for an ogive cylinder, with a plane perpen- 
dicular to the axis of the ogive cylinder acting as a two-dimensional patch boundary 
separating two three-dimensional regions agree well with experimental data. The 
factored, implicit patched-grid approach was used to develop the three-dimensional 
code . v 

The flux balance in three dimensions is performed on a surface instead of a 
line (CD in fig. 8). In order to calculate the interpolation coefficients Nj 
the overlap areas of the control volume sides that come together along the flux 
balance surface need to be calculated. In the case of planar patch boundaries, this 
amounts to finding the area of overlap of two triangles given the locations of their 
vertices. This can be a time consuming process depending on the method used. 
Investigations to find the simplest way of calculating the area of overlap of two 
general triangles are currently in progress. The interpolation required to enforce 
the continuity of dependent variables across the patch boundary also needs to be 
two-dimensional in nature to perform three-dimensional patched-grid calculations. 

The trial calculation mentioned above uses linear interpolation over triangles. 
Details of three-dimensional patched-grid calculations will be discussed in a forth- 
coming article. 


Patched-Grids and Navier-Stokes Calculations 

The Navier-Stokes equations can be written in the conservation-law form (just 
as the Euler equations were in eq. (1)), with the use of viscous flux vectors. 

These viscous flux vectors can be balanced across the flux balance line in a manner 
identical to the inviscid flux balance discussed in the last section. The viscous 
fluxes depend on both the dependent variables and their spatial derivatives. These 
fluxes can be evaluated using the dependent variables and their derivatives on one 
side of the patch boundary, and used to update the dependent variables on ohe other 
side of the patch boundary. 

However, in many situations a much simpler technique may suffice. Most viscous 
calculations in the area of applied aerodynamics are performed with the thin-layer 
Navier-Stokes equations. The viscous terms in the direction parallel to the body 
surface are assumed to be negligible in comparison with the viscous terms in the 
direction normal to the body surface to obtain the thin-layer equations. Even in 
cases where this assumption does not hold true, the computer speeds and memory 
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available currently preclude the use of very fine grids in more than one spatial 
direction. Coarse grids will not resolve viscous terms and, hence, retaining vis- 
cous terms in a direction in which the grid is coarse does not serve any purpose 
(the magnitude of the viscous terms calculated on such grids will be negligibly 
small compared to the other terms). Most practical applications will require the 
use of patch boundaries that intersect the surface boundary and not patch boundaries 
that closely follow the surface (or are parallel to the surface). The use of the 
thin-layer equations in such situations does not require a flux balance of the 
viscous terms. However, the patch-boundary condition can be made fully conserva- 
tive, with very few modifications, for the complete Navier-Stokes equations if the 
need arises. 


Higher-Order Accurate Patch-Boundary Conditions 

The linear interpolation used to enforce continuity of the dependent variables 
across the patch boundary and the piecewise-constant variation of the numerical 
fluxes assumed for the initial development of the patch-boundary condition lead to 
first-order accuracy at the boundary (even if the interior scheme is second-order 
accurate in nature). Further, a sudden change in spacing (spacing normal to the 
patch boundary) across the patch boundary could lead to additional deterioration of 
the solution accuracy in the vicinity of the patch boundary. To overcome these 
problems, it is necessary to use higher-order accurate interpolation schemes and a 
piecewise linear (or more accurate) variation of the numerical fluxes. In addition, 
the accuracy with which numerical fluxes are calculated must be made independent of 
sudden jumps in spacing (in the direction normal to the patch boundary). Higher- 
order accurate interpolation schemes must be used with care because they may result 
in unstable patch-boundary procedures. Higher-order accurate patch-boundary proce- 
dures are currently being investigated and will be discussed in a future article. 


RESULTS 


Results are presented in this section for inviscid supersonic flow past an 80° 
segment of a cylinder, subsonic flow past a full cylinder, blast-wave diffraction by 
a ramp, vortex motion through a patch boundary and time-dependent flow past airfoils 
that move relative to each other. The unsteady Euler equations are integrated in 
time in all these cases. The results presented in this section were obtained with 
the explicit, factored implicit and implicit-relaxation schemes and give a flavor of 
the quality of solutions possible with each scheme. The case of supersonic flow 
over a cylinder demonstrates the shock-capturing quality of the various schemes and 
the convergence rates possible with each one. The blast-wave simulation shows the 
applicability of patched grids in situations requiring selective grid refinement and 
the vortex calculation shows the time-accuracy of the patch-grid schemes and the 
distortion-free motion of vortices through patch boundaries. The subsonic cylinder 
case demonstrates the time-accuracy of the implicit-relaxation patch-boundary scheme 
and the feasibility of performing calculations on patches that move relative to each 
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other. The double-airfoil calculation shows the applicability of the moving-patch 
technology to rotor-stator interaction calculations. 


Cylinder in a Supersonic Free Stream 

The first test case is that of a cylinder in a supersonic free stream (M m = 2) 
with the associated bow shock. Figure 10 shows the grid used for the calculation. 
The region of interest was divided into two patches separated by the patch boundary 
AB. The discontinuity of the constant 5 grid lines at the patch boundary is 
evident. The values of the dependent variables at all the grid points were set 
equal to their free stream values initially. The equations of motion together with 
the various boundary conditions were integrated (including the patch-boundary condi- 
tions) until the solution converged to its steady-state value. 

Figure 11 shows the pressure contours obtained at convergence with the 
explicit, first-order accurate Osher scheme and the explicit patch-boundary 
scheme. The contour lines can be seen to be continuous across the patch boundary 
(AB). In fact the contour lines seem almost slope continuous across AB. The square 
symbols in this figure (and the following figures pertaining to the cylinder) repre- 
sent the shock position predicted by another numerical approach in reference 26. 

The captured shock is a little to the left of the predicted shock. This discrepancy 
is characteristic of the first-order accurate Osher scheme and will disappear with 
the use of a second-order accurate integration scheme. 

The bow shock associated with the supersonic flow over a cylinder first appears 
at the surface of the cylinder and then moves outward to its final converged posi- 
tion (for the initial conditions chosen). Figures 12-15 depict pressure contours 
showing the progression of the shock through the grid system. The results in fig- 
ures 12-15 were obtained with the first-order accurate, implicit, relaxation scheme, 
and the corresponding implicit patched-grid scheme. The pointwise relaxation tech- 
nique was employed for this calculation. Because of the large transients that occur 
during the first few steps of the calculation, the CFL number was initially 
restricted to 10.0 (the first 10 steps) and thereafter increased to 500.0. The use 
of CFL numbers larger than 500.0 did not alter the convergence rate. 

Figure 12 shows the pressure contours after 20 integration steps. The shock at 
this instant in time is quite distant from the patch boundary. Figure 13 shows the 
pressure contours after 35 steps, and figure 14 the contours after 60 steps. The 
smooth transition of the shock from patch 1 to patch 2, even at a CFL number of 
500.0, is evident. Figure 15 presents the pressure contours obtained at convergence 
(after 280 steps). The continuity of the contour lines across the patch boundary in 
figures 14 and 15 demonstrate the quality of solutions possible with the present 
patch-boundary scheme. 

Figure 16 displays the convergence history for the explicit first-order scheme 
and the implicit relaxation first-order scheme. The explicit scheme required 
approximately 2700 steps to converge. The convergence criterion chosen for this and 
the following cylinder calculations was 
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(based on a freestream density value of 1.0). The relaxation scheme with only one 
iteration per step required 280 steps to converge; that is, the use of the implicit 
relaxation scheme increased the convergence rate by a factor of 9.65. Since the 
implicit scheme required 1.22 times as much computing time per step as the explicit 
scheme, the actual computing cost was reduced by a factor of 7.91 with the use of 
the relaxation scheme. It should be remembered that the extra programming required 
to implement pointwise relaxation is minimal (typically, less than 100 lines of 
Fortran). Figure 16 also shows the convergence rate for the implicit relaxation 
scheme when two_ iterations are used at each time step. The scheme in this case 
required only 157 steps to converge. However, since it requires approximately twice 
as much computing time per step (compared to the relaxation scheme with one itera- 
tion per step) the total computing cost was almost the same. Figure 17 shows the 
variation of the magnitude of the maximum residue (R) in the continuity equation as 
a function of the number of integration steps 
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where the superscript c represents the element of the vector corresponding to the 
continuity equation. The relaxation scheme with one iteration per step reduces R 
by 9 orders of magnitude in about 1000 steps and the scheme with two iterations per 
step required about 500 steps to decrease R by the same amount. Table 1 summar- 
izes the convergence rates and computing times for the variants of the first-order 
accurate Osher scheme mentioned above. All computations were performed on a CRAY- 
XMP. The computer programs used in the present study were vectorized for use on 
this machine. Additional information regarding computing times and convergence 
rates for other relaxation techniques such as the line-relaxation approach can be 
found in reference 9. 

First-order accurate schemes are insufficient to provide accurate results for a 
general class of problems. In order to obtain reliable results it is necessary to 
resort to second-order accurate integration schemes. The first-order results pre- 
sented in this study were included merely to demonstrate the increase in convergence 
rate (compared to that obtained with the explicit patched-grid scheme). Figure 18 
shows the pressure contours obtained at convergence with the implicit, relaxation 
second-order accurate Osher scheme for the cylinder problem. This scheme uses flux 
limiters to achieve the TVD property in each dimension. The shock position can be 
seen to be predicted more accurately than in the previous case. The contours once 
again transition smoothly across the patch boundary despite the discontinuity of the 
grid lines. 

Figures 19-22 are very similar to figures 12-15 in that they show the progres- 
sion of the shock through the grid system to its final converged position. However, 
the results shown in figures 19-22 were obtained with the implicit, factored, 
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first-order accurate Osher scheme. Once again, because of the large transients that 
occur in the initial stages of the calculation, the CFL number was restricted to 5.0 
in the first 10 steps, and was then held constant at a value of 40.0 for the rest of 
the calculation. Unlike the implicit relaxation scheme, the implicit factored 
scheme became unstable for CFL numbers larger than 40.0 for this case. Figures 19, 
20, and 21 show the pressure contours obtained after 13, 19, and 24 integration 
steps, respectively. The shock is seen to be passing through the patch boundary in 
figure 20, and in figure 21 the shock has passed through the patch boundary. 

Figure 22 shows the pressure contours at convergence. 

Unlike the shock transition, from patch 1 to patch 2 obtained with the 
implicit, relaxation scheme, the shock transition obtained with the implicit- 
factored scheme is not distortion free. This can be clearly seen in figures 20 
and 21. The oscillations near the shock are seen even before the shock encounters 
the patch boundary. These oscillations grow larger as the shock passes through the 
boundary and eventually vanish at convergence (fig. 22) as expected. This phenome- 
non is because of the inconsistency between the linearization of the fluxes at the 
patch boundary and the fluxes in the interior of each patch employed in the fac- 
tored, implicit patched-grid scheme. Such an inconsistency is absent in the case of 
implicit, relaxation patched-grid schemes where the pointwise relaxation technique 
or the line-relaxation technique (with the lines chosen to belong to the same family 
as the patch boundary) are employed. Since the distortion is because of 
linearization error alone, it vanishes when the solution converges or if the itera- 
tion procedure is carried to convergence at each time step. To demonstrate the 
restoration of the properties of the fully implicit scheme (when the iterations 
converge), the following test was performed. The solution obtained after 22 inte- 
gration steps (from the previous calculation) was used as the starting solution and 
was integrated through two time steps. A CFL number of 40.0 with 5 iterations per 
step was used. The residual was reduced by approximately an order of magnitude at 
the fifth iteration. The pressure contours obtained as a result of this calculation 
are presented in figure 23. Clearly, the oscillations of figure 21 are absent in 

figure 23 (figs. 21 and 23 can be directly compared since they show pressure con- 

tours at the same time step). The monotonicity property is restored, because at 
convergence, the difference equations corresponding to the fully implicit scheme are 
satisfied. 

Figure 24 shows the convergence history for the factored, implicit, first-order 
calculation (with the implicit patch-boundary condition) and for the explicit first- 
order calculation (with the explicit patch-boundary condition). The convergence 
criterion was chosen to be the same as before. The factored, implicit sch.me 
required 94 steps to converge, whereas the explicit scheme required 2700 steps to 
converge; that is, the use of the factored, implicit scheme increased the overall 
convergence rate by a factor of 28.7. However, the factored, implicit scheme with 
two iterations per time step (the results shown in figs. 19-22 were obtained with 

two iterations per time step) requires 3.35 times as much computing time per time 

step as the explicit scheme. Hence, the actual computing cost was reduced by a 
factor of 8.6 with the use of the factored, implicit scheme. Figure 25 shows the 
variation of the magnitude of the maximum residual in the continuity equation (R) as 
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a function of the number of integration steps. The factored, implicit scheme with 
two iterations per step reduced R by 9 orders of magnitude in about 250 steps. 

The explicit scheme requires about 3000 steps to reduce R by two orders of 
magnitude. Table 2 summarizes the computation times and convergence rates for the 
factored implicit scheme. 

Figure 26 displays the pressure contours obtained at convergence with the 
factored, implicit, second-order accurate Osher scheme. Unlike the implicit, relax- 
ation scheme, the factored implicit scheme does not require the TVD property to be 
satisfied in each spatial direction. Hence, the fluxes used for the calculation 
shown in figure 26 were not postprocessed using a flux limiter. For this reason the 
pressure contours in the vicinity of the shock are oscillatory. However, the oscil- 
lations are confined to a small region-near the shock.. 

A general-purpose, patched-grid Euler code should have the capability of han- 
dling as many patches as necessary to cover the region of interest. In order to 
demonstrate the generality of the present patched-grid scheme and its applicability 
in a general-purpose, patched-grid Euler code, the region of interest for the cylin- 
der was divided into four patches instead of two patches as in the previous case. 

The patches and the grids for each patch are shown in figure 27. A wide variation 
in cell shapes and sizes can be seen across each of the three patch-boundaries. 
Figure 28 shows the pressure contours obtained at convergence with the factored, 
implicit, first-order accurate Osher scheme. Small changes in contour line slopes 
can be seen at the patch boundaries. This is because of the sudden change in accu- 
racy caused by the abrupt transition in mesh clustering and will be less noticeable 
for higher-order accurate schemes. Figure 29 shows the pressure contours obtained 
with the factored, implicit, second-order accurate Osher scheme. Once again, the 
contours are continuous across the patch boundaries (and slope continuous unlike the 
contours in fig. 28). The smooth transition of the shock from patch 2 to patch 4 
emphasizes the quality of solutions possible with the patched-grid schemes presented 
here. 


Blast-Wave Diffraction by a Ramp 

As stated earlier, one of the main advantages of being able to perform patched- 
grid calculations is that one can selectively refine the grid in certain areas of 
the flow region without having to maintain grid-line continuity across patch bound- 
aries. The problem of blast-wave diffraction by a ramp with its attendant complex 
shock patterns is a typical example of a problem requiring selective refinement of 
the grid. The grid used for the calculation presented here is shown in figure 30. 
The calculation was performed with an incident shock Mach number of 7.1, a ramp 
angle of 49° and a ratio of specific heats (y) of 1.55. A double Mach reflection 
occurs for this choice of flow parameters. The discontinuities for this configura- 
tion include the incident shock, a kinked reflected shock, two Mach stems, and two 
slip lines. The region in which the triple point and Mach stems occur (the area 
covered by patch 2) requires a very fine grid in order to resolve all the flow 
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features. A very coarse grid is sufficient for the area covered by patch 3 because 
the independent variables are constant in this region. 

The solution to the blast-wave diffraction problem is self similar in time and, 
hence, can be made to have an asymptotic limit in computational space with the help 
of a similarity transformation. Details of the similarity transformation can be 
found in reference 27. Figure 31 shows the isobars obtained at convergence with the 
explicit, first-order accurate Osher scheme. The kinked reflected shock (ABC), the 
incident shock (CD), and the first Mach stem (CE) are seen to be captured very 
sharply. In addition, the second Mach stem (much weaker than the first) can be seen 
emanating from the point B. The reflected shock passes through the patch boundary 
FG without any distortion; it merely becomes thinner in patch 2 because of the fine 
grid used in this patch. Figure 32 shows the density contours obtained at conver- 
gence. The additional feature that this figure brings to light is the first slip 
line emanating from the triple point C. The second slip line, which is supposed to 
emerge from point B, is much weaker than all the other discontinuities, and hence, 
cannot be seen in this figure. 


Vortex Motion Through a Patch Boundary 

This calculation consisted of a Lamb-type analytical vortex moving through a 
patch boundary. It is possible to affect vortex motion either by superimposing a 
moving free stream condition (in which case the vortex is convected along with the 
fluid at the free stream velocity), or by keeping the vortex stationary and moving 
the grid. The two approaches yield identical results. In this study the vortex was 
initialized using the procedure given in reference 28 and then the grid was moved at 
a constant speed in a direction opposite to the direction in which vortex motion was 
desired. Figure 33 shows the two-patch grid used for the calculation. Only the 
central portion of the flow field is presented in figures 33-38 since the essential 
features of the vortex are contained in this region. The discontinuous nature of 
the grid lines at the patch boundary is evident.- The calculation was performed with 
the factored, implicit, second-order accurate Osher scheme and three iterations per 
step. Second-order accuracy in time was achieved by using a three-point backward 
difference for the time derivative on the right-hand side of equation (23). 

Figure 34 shows pressure contours at- initialization. The solid core at the , 
center of the constant pressure circles in this figure and the following figures 
that show pressure contours, is the analytical center of the vortex. Figure 35 
shows pressure contours after 85 integration steps. The slope continuity of the 
contour lines across the patch boundary is clearly seen. The center of the circles 
coincides with the analytical center of the vortex (this demonstrates the time- 
accuracy of the patched-grid scheme). Figure 36 presents pressure contours obtained 
after 1 80 integration steps. The vortex has moved entirely into patch 2. The 
contours are circular and are not distorted. However, because of asymmetric trunca- 
tion errors the center of the circles is slightly above the analytical center of the 
vortex. Figure 37 shows a continuous grid on which the calculation was performed 
once again. The results of this computation are shown in figure 38 (pressure 


31 



contours after approximately 140 integration steps). Once again, a slight upward 
movement can be observed, demonstrating that this motion is not due to an inaccurate 
patch-boundary condition. 


Cylinder in a Subsonic Free Stream (Moving Patches) 

This calculation was performed for a complete cylinder in subsonic flow. The 
free stream Mach number for the calculation was 0.35. Figure 39 shows the two grids 
used to perform the calculation. All the grid points were initialized to their free 
stream values, and the governing equations together with the patch-boundary condi- 
tion were integrated to convergence (both the grids were stationary for this phase 
of the calculation). The integration scheme used was the factored, implicit, Beam- 
Warming scheme in a patched-grid setting (ref. 7). Second-order accuracy in time 
was achieved as in the previous case. The outer grid was then made to rotate at a 
constant angular speed. An integration scheme without any truncation error would 
result in the dependent variables being stationary in physical space, and conse- 
quently, contour levels of pressure, density, etc., would also remain stationary 
(although in computational space, the dependent variables at a grid point of the 
outer grid would be constantly changing because of the changing physical location of 
the grid point). 

Figure 40 shows the pressure contours obtained at each time step, as the outer 
grid performed one rotation, plotted on the same figure. The number of integration 
steps used per rotation was 260. Five iterations were performed at each step in 
order to reduce the magnitude of the residue by approximately an order of magni- 
tude. For an integration scheme of infinite accuracy, the thickness of the bands 
would be zero. Since the integration scheme used is only second-order accurate in 
space and time, the thickness of the contour bands is finite. The thickness of the 
bands gives a qualitative measure of the time accuracy of the integration scheme 
coupled with the patched-grid scheme. Further details of this calculation can be 
found in reference 29. Reference 29 also includes a demonstration calculation which 
shows that the finite thickness of the contour bands is because of the accuracy of 
the integration scheme and not because of any inadequacy of the patch-boundary 
scheme. 


Time-Dependent Double-Airfoil Calculation 

One of the important applications of the patched-grid approach is the treatment 
of regions that contain bodies which move relative to other bodies such as rotor- 
stator combinations in turbomachinery. The feasibility of performing such calcula- 
tions has already been demonstrated in the previous calculation simulating flow over 
a -full cylinder. In this subsection, results obtained with the unsteady Euler 
equations for a simple rotor-stator configuration in which both the rotor and stator 
are circular arc airfoils, are presented. The axial gap between them is 20 % of the 
chord length. The free stream Mach number is 1.5. The integration scheme used is 
the implicit relaxation second-order accurate Osher scheme. Second-order accuracy 
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in time is obtained as in the previous case. Flux limiters have been used to make 
the scheme TVD in each spatial dimension. 

Figure 41 shows the two-patch grid used to discretize the region of interest. 
The grid in patch 2 is stationary (as is the aft airfoil), and the grid in patch 1 
is fixed to the first airfoil which is moving vertically downward. Although the 
grid lines are continuous at the patch boundary in figure 41, a discontinuity in 
these grid lines will develop as the first airfoil and patch 1 move downward. The 
patch-boundary points of patch 1 will slip past the patch-boundary points of 
patch 2. 

Periodic boundary conditions are imposed on the upper and lower boundaries of 
both patches. Free stream conditions are imposed on the left boundary of patch 1, 
and supersonic exit boundary conditions are imposed on the right boundary of 
patch 2. The implicit, patch-boundary condition is used at the boundary separating 
the two patches. Details of the inviscid surface boundary condition used can be 
found in reference 10. 

The calculation was initially performed with both airfoils stationary. After 
50 integration steps, the first airfoil was given a downward velocity (the magnitude 
of the velocity corresponding to a Mach number of 0.1 with respect to free stream 
conditions), and the calculation was continued until the solution became periodic in 
time. The calculation was performed at a CFL number of approximately 2.0. At this 
CFL number, 250 integration steps were required for each cycle (one cycle 
corresponds to the motion of the upper boundary of patch 1 from its current position 
to the position occupied currently by the lower boundary of patch 1). Approximately 
three cycles were required to eliminate the initial transients and to establish a 
solution that was periodic in time. 

Figures 42-47 show pressure contours at various positions of the forward air- 
foil (with respect to the aft airfoil) as it moves downward. These contours were 
obtained after the initial transients had subsided. Although the calculation was 
performed with only two airfoils, for the sake of clarity figures 42-47 depict four 
airfoils (information regarding the additional airfoils was obtained from the peri- 
odicity condition). Figure 42 shows contours at t = 0.0 (one cycle corresponds 
to t = 1.0). The downward motion of the forward airfoil results in an effective 
angle of attack which, in turn, results in an attached oblique shock on the lower 
side, and a weak expansion fan on the upper side at the leading edge of the first 
airfoil. The interaction of this shock with the adjacent forward airfoil is clearly 
seen. A second, weak attached shock is also evident at the trailing edge (lower 
side) of the forward airfoil. At the position t = 0.0, the shock associated with 
the second airfoil is detached. It is seen impinging on the surface of the adjacent 
airfoil. 

The interaction of the trailing edge shock of the first airfoil and the leading 
edge shock of the second airfoil is also clear from figure 42. This interaction 
area moves downward as the first airfoil moves downward. This, in turn, results in 
the leading edge shock of the second airfoil attaching and detaching periodically 
from the leading edge. Figures 43-47 depict this attachment/detachment process. In 
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figure 43 the shock is beginning to reattach. In figure 44 it is a weak attached 
shock, and in figure 45 it is a strong attached shock (beginning to detach again). 

It is completely detached in figure 46 and finally, in figure 47 the contours are 
identical to those in figure 42 thus demonstrating the accuracy of the present 
technique in calculating periodic flows. 

An important feature in figures 42-47 is that the contours are continuous, even 
slope continuous, across the patch boundary along which grid points from the two 
grids are slipping past each other. This high degree of continuity is because of 
the conservative nature of the patched-grid scheme, and the manner in which continu- 
ity of dependent variables is enforced across the patch boundary. (Another inter- 
esting feature is that the captured shocks are almost oscillation-free. The absence 
of large oscillations is because of the TVD nature of the integration scheme.) 

Figure 48 shows the pressure variation with time at midchord on the lower 
surface of the aft airfoil. This pressure variation corresponds to the fourth and 
fifth cycles. Clearly, the pressure is periodic in time, thus demonstrating the 
capability of the integration, and patch-boundary schemes in simulating periodic 
time-dependent flow. Figure 49 shows the surface pressure variation at midchord on 
the upper surface of the aft airfoil. The behavior seen in figure 49 is similar to 
that seen in figure 48 except for a phase shift and a difference in the mean value 
of the pressure. 

One aspect of rotor-stator configurations that is not represented in the pres- 
ent results is the effect of the aft airfoil on the forward airfoil (the supersonic 
nature of the flow does not permit such an interaction). However, the patch- 
boundary conditions presented earlier will accurately simulate such an interaction, 
if it were present. The calculation of reference 23 (presented in the next section) 
with purely subsonic flow in the region of interest shows the capability of the 
patched-grid scheme in accurately predicting such an interaction. 


AM APPLICATION TO ROTOR-STATOR INTERACTION 


The aerodynamic processes associated with the flow of fluid through turbo- 
machines pose one of the most challenging problems for the computational fluid 
dynamicist. The unsteady nature of the flow, the complex geometries involved, the 
motion of some parts of the system relative to others, and the periodic transition 
of the flow from laminar to turbulent are some of the factors that contribute to the 
complexity of the problem. A clear understanding of these types of flows is essen- 
tial for the optimization of the performance of turbomachinery. This section pre- 
sents an analysis of two-dimensional flow past the rotor-stator configuration of an 
axial-turbine (ref. 23), using the patched-grid methodology described earlier. 

The two-dimensional analysis of stator airfoils in isolation or rotor airfoils 
in isolation is a relatively straightforward task. Such an analysis is valid when 
the two rows of blades are set far enough apart that the interaction effects are 
minimal. However, the desire to minimize engine length requires the stator and 
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rotor rows to be closely spaced. Clearly, the interaction effects will become more 
important as the axial gap between the rows is reduced. In fact, the flow becomes 
periodically unsteady for small values of the axial gap. The experimental results 
of reference 30 show that the pressure fluctuation (difference between the minimum 
and maximum pressure values) near the leading edge of the rotor airfoil can be as 
much as 12 % of the exit dynamic pressure when the axial gap is reduced to 15$ of the 
chord length (for the operating conditions and geometry chosen). Hence, it is 
important that the rotor and stator airfoils be treated as a single system when the 
interaction effects become predominant. A computational tool which provides the 
design engineer with the necessary aerodynamic data can be used to great advantage 
in redesigning rotor and stator airfoils to enhance performance. Such a tool would 
have to accurately simulate the unsteady flow associated with rotor-stator config- 
urations which exhibit a strong interaction. 

The multiple-grid approach can be used to great advantage in simulating the 
flow associated with rotor-stator configurations. The motion of the rotor airfoil 
with respect to the stator airfoil makes it impractical to wrap a single grid around 
both airfoils. It is much simpler to contain the stationary airfoils in one set of 
grids (grids that are stationary) and the moving airfoils in another set of grids 
(grids that move along with the airfoils) and to have the two sets of grids communi- 
cate with each other using patched-grid technology. The ability to accurately 
transfer information from patch to patch, in such a situation, has already been 
demonstrated in the previous section. 

One aspect of rotor-stator configurations that is not represented in the rotor- 
stator results of the previous section is the effect of the aft airfoil on the 
forward airfoil (the supersonic nature of the flow does not permit such an interac- 
tion). However, the patch-boundary conditions were implemented such that an inter- 
action, if present, would be properly accounted for. The calculation presented in 
this section, with the blade geometry of reference 30 and purely subsonic flow in 
the region of interest, tests the capability of the patched-grid scheme in accu- 
rately simulating such an interaction. 

Although the relative merits of overlaid and patched grids have been discussed 
in the introductory section, the procedures used to transfer information between 
overlaid grids have not been presented. The main advantage of overlaid grids is 
that, unlike patched grids, the outer boundaries of the individual grids do not have 
to conform to each other. This aspect of overlaid grids can be very useful when 
dissimilar geometries are to be coupled; for example, when a circular body is 
required to fit into a rectangular flow region as in figures 1(a) and 1(b) (the 
overlaid approach requires fewer zones). Reference 31 deals with the use of over- 
laid grids in solving the Euler equations. Encouraging results are presented for 
several example problems. However, the issues of time-accuracy and conservation are 
not addressed in reference 31. This section includes the necessary extensions to 
the overlaid grid technique of reference 31 in order to make it time-accurate. 

In this section, both the patched and overlaid grid techniques are used to 
simulate the flow associated with the rotor-stator configuration of reference 30. 

The thin-layer, time-dependent, Navier-Stokes equations in two spatial dimensions 
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are solved using a variant of the Beam-Warming scheme (ref. 5) that is developed in 
reference 8. Comparisons with the experimental data of reference 30 are also 
made. The comparisons include time-averaged values of surface pressure, pressure 
amplitudes and phase relationships. 


Integration Method 

The integration method used is the factored, implicit algorithm developed 
earlier with modifications to include the viscous terms. To describe the scheme we 
consider the unsteady, Navier-Stokes equations in two dimensions 


Q. + E + F = R +S (49) 

t x y x y v ' 

where the vectors Q, E, and F are given in equation (2). The vectors R and S 
are given by 
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Making the independent variable transformation given in equation (4) and the 
thin layer approximation (ref. 32), we obtain 


Q + E + E = Re _1 S 
t 5 n n 


(52) 


where Q, E, and F are given in equation (6) and the vector S is given by 
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Equations (52) and (53) assume that the body surface is a constant n line (in 
incorporating the thin-layer approximation). 

The factored, implicit algorithm for the unsteady, thin-layer, Navier-Stokes 
equations is given by 
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(54) 


A* = ( BE/aQ) 1 
B ± = (aF/aQ) 4 
M = aS/aQ 


and A, V, and 5 are forward, backward, and central difference operators, respec- 
tively. The quantities E j + i /2 E j k+1/2’ and S j k-tl/2 are numerical fluxes 
consistent with the physical fluxes £, F, and S. When second-order accuracy in 

time is required, the term (Q^ , - Q 1 . 1 , ) on the right-hand side of equation (54) 

J > k J , k 

must be replaced by (1.5Q^ . - 2.0Q a . + O.SQ 1 ?"*!) in addition to iterating to con- 

J ) K J ? K J ) 

vergence. Typically, three to four iterations per time step are sufficient to 
reduce the residual by an order of magnitude or more. 
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Different choices for the numerical fluxes lead to different integration 
schemes. The central difference scheme used in this study can be obtained with the 
numerical fluxes defined as 


1 rS 


where 


E j+1/2,k = 2 [E(Q j ,k'^j+1/2,k ) + E(Q j+1,k’ 5 j+1/2,k )] ) 


E j , k+1/2 = 2 [E(Q j,k’ n j,k+1/2 ) + ?(Q j,k+1 ,n j,k+1/2 )] 


S j , k+1/2 = §[Q j,k+1/2' (Q n ) j,k+1/2’ T1 j,k + 1/2 ] 


Q j , k+1/2 = 2 (Q j,k + Q j,k+1 } 


(55) 


(Q n ) j, k+1/2 = Q j,k+1 " Q j,k 


The smoothing terms necessary to stabilize the calculation were included in the 
numerical fluxes by redefining them to be 


E j+1/2,k = 2 ^ E(Q j,k ,? j+1/2,k ) 


E(Q j+1 ,k’^j+1/2,k ) 


+ 2 [hE(Q j _ 1)k> Q j>k ,C j+1/2>k )| - 2 |A g (Qj jk ,Qj +1fk ,5j +1/2>k ) I 


l AE(Q j+1,k ,Q j+2,k’ 5 j+1/2,k ) l ] ^ 


(56) 


where 


I AE(Q j ,k’ Q j+1 ,k’^j+1/2,k ) I = {Q J+1/2,k^j+1/2,k ) l (Q j+1,k " 5 j,k } (57) 


The numerical flux Fj k+ ']/ 2 i- s redefined in a similar manner. The advantage in 

including the smoothing term in the definition of the numerical flux is that it can 
also be made conservative across patch boundaries. The term e in equation (56) 
determines the amount of smoothing used in the calculation, e = 1 corresponding to 
the amount of smoothing in a fully upwind scheme. Additional details of this 
smoothing parameter can be found in reference 7. 
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Grid System for the Rotor-Stator Configuration 


The airfoil geometry of reference 30 is much more complicated than the simple 
circular-arc airfoils used in the previous section. Figure 50 shows the shape of 
the rotor and stator airfoils, and their orientation relative to each other and the 
inlet flow direction. The stator and rotor rows consisted of 22 and 28 airfoils, 
respectively. A Navier-Stokes calculation with a total of 50 airfoils would be 
extremely time consuming. To overcome this difficulty, the rotor airfoils were 
enlarged by a factor of (28/22) keeping the pitch to chord ratio the same. It was 
then assumed (for the calculation) that there were only 22 airfoils in the rotor 
row. This assumption makes it possible to perform a calculation with only one rotor 
airfoil and one stator airfoil (a periodicity boundary condition is used to simulate 
the presence of the rest of the airfoils). The original stator airfoil and the 
modified rotor airfoil are shown in figure 51. The axial gap between them was 
chosen to be \ 5 % of the average chord length. 

This rotor-stator configuration and the associated flow region can be discre- 
tized using only patched grids. However, the number of patches required to solve 
the problem accurately would be twice as many as required when both patched and 
overlaid grids are used in conjunction with each other. Hence, a combination of 
patched and overlaid grids is used in this calculation. 

The grid system chosen for the calculation consists of four different zones. 

The first zone contains the stator airfoil and is discretized with an "0" grid. The 
second zone contains the rotor airfoil and is also discretized with an "0" grid. 
These two grids do not overlap one another. Grids 1 and 2 were generated using an 
elliptic, grid generator of the type developed in reference 33, and they are shown in 
figure 52. The set of grid lines that intersect the airfoil surfaces are orthogonal 
to these surfaces. Although the actual grids used for the calculation are very 
dense near the airfoil surfaces (to resolve the viscous effects), for the purpose of 
clarity figure 52 shows grids in which the points are equispaced in the direction 
normal to the airfoil surfaces. 

Grids 3 and 4 are shown in figure 53. These grids are generated algebrai- 
cally. Grid 3 contains grid 1 and, consequently, the stator airfoil. Grids 1 and 3 
overlap each other and are also stationary with respect to each other. In fact, the 
inner boundary of grid 3 is contained within grid 1 and the outer boundary of grid 1 
is contained within grid 3. This positioning of the boundaries is necessary for the 
transfer of information between grids 1 and 3. The relationship between grids 2 
and 4 is similar to that between grids 1 and 3. Additionally, grids 3 and 4 abut 
each other along the patch-boundary ABCD (fig. 53). These two grids (3 and 4) slip 
past each other as the rotor airfoil moves downward. Grids 3 and 4 constitute a 
patched-grid system. It is advantageous to use a patch boundary (as opposed to an 
area of overlay) where one system of grids moves relative to another because both 
time-accuracy and conservation can be more easily controlled in patched-grid 
calculations . 

An interesting feature of grids 3 and 4 as seen in figure 53 is that they do 
not align with each other. The segment AB of grid 4 does not seem to align with any 
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part of the patch-boundary of grid 3, and similarly, the segment CD of grid 3 does 
not seem to align with any part of the patch boundary of grid 4. However, the 
periodicity boundary condition can be used to solve this problem; the result being 
that the segment AB is matched with the segment CD. Figure 54 shows all the grids 
used in the calculation. The relative positions of the boundaries of each grid are 
clearly seen. 


Boundary Conditions 

The use of_ multiple grids in simulating the flow over the rotor-stator config- 
uration shown in figure 51 results in several computational boundaries (fig. 54). 

The boundary conditions used at each of these boundaries is briefly outlined below. 

Natural boundaries - The inner boundaries of the two "0" grids correspond to the 
airfoil surfaces and, hence, the "no slip" condition and adiabatic wall conditions 
(or wall temperature) are imposed at these boundaries. It should be noted that in 
the case of the rotor airfoil "no slip" does not imply zero velocity at the surface 
of the airfoil; instead, it means that the fluid velocity at the rotor surface is 
equal to the rotor speed. In addition to the "no slip" condition, the derivative of 
pressure normal to the wall surface is set to zero. The pressure derivative condi- 
tion, the adiabatic wall condition and the equation of state together yield 


3p 

an 


0 


ae _ u 3pu v apv 

an " p an p an 


where n is the direction normal to the airfoil surface. These boundary conditions 
are implemented in an implicit manner by using the following equation instead of 
equation (54) to update the grid points on the airfoil surfaces 


where 


c(C! - ) + D(QP + ^ - QP ,) = 0 


\5,1 "5,1 


J,2 "j,^ 


c 



0 0 
1 0 
0 1 
a e 




(58) 
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and 
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8 



u wall 

p wall 

V wall 

p wall 


Equation (58) is an implicit, spatially first-order accurate implementation of the 
"no slip" and adiabatic wall conditions (first-order accurate because the zero 
normal derivative condition is implemented using a two-point forward difference). A 
second-order accurate, three-point forward difference corrector step is also imple- 
mented after each time step. It should be noted that equation (58) requires the 
grid to be orthogonal at the airfoil surfaces and the Jacobians of the transforma- 
tion J. * and J , o to be independent of t. 

J 9 1 J t 

The left boundary of grid 3 is a subsonic inlet boundary. Three quantities 
need to be specified at this boundary. Those chosen for this study are (refs. 34 
and 35) the Riemann invariants 


R 


1 


u + 




(59) 


and the inlet flow angle which in this case is equivalent to 


v inlet 


( 60 ) 


The fourth quantity (which is necessary to update the points on this boundary) is 
also a Riemann invariant 


R„ = u - — (61) 

d. y - 1 

and is extrapolated from the interior of grid 3. The manner in which these boundary 
conditions can be implemented implicitly is described in reference 36. 

The calculation assumes that there are an infinite number of rotor and stator 
airfoils in the positive and negative y directions in figure 51. Hence, a, simple 
periodicity boundary condition is imposed on the upper and lower surfaces of grids 3 
and 4. The implicit implementation of this boundary condition is straightforward 
and will not be discussed here. The right boundary of grid 4 is a subsonic exit 
boundary. A simple extrapolation procedure is used at this boundary. This boundary 
condition is implemented implicitly by using the equation 


41 



( 62 ) 


_ 3 P x _ iimax-1_ ,k ,*p+1 _ *P ) _ 0 

,k jmax,k J JmaX)k jmax-1,k jmax-1,k' 

instead of equation (54) to update the grid points. While the use of equation (62) 
results in a stable calculation, it may be necessary to specify one condition (such 
as R 2 or the exit static pressure) to obtain accurate results in the vicinity of 
the exit boundary. The importance of such a modification to the exit boundary 
condition is currently being investigated. 

Overlap boundaries - The overlap boundary conditions used at the inner bound- 
aries of grids 3 and 4 and the outer boundaries of grids 1 and 2 are extremely 
simple. They are given by 

(Q P+1 - Q P ) 0>b> = 0 (63) 


(Q p+1 

jmax 


where the subscript o.b. refers to the points on the overlap boundary. This is 
followed by an explicit, corrective, interpolation procedure at the end of each 
iteration wherein the values of Q p+1 along the outer boundaries of grids 1 and 2 
are interpolated from the interior grid points of grids 3 and 4, respectively, and 
the values of Q p+1 along the inner boundaries of grids 3 and 4 are interpolated 
from the interior grid points of grids 1 and 2, respectively (the values of Q p+1 
obtained from the implicit integration procedure along the overlap boundaries are 
discarded). The results in this study were obtained using a linear interpolation 
over triangles. 

The boundary condition given by equation (63) is not the same as that given by 
the following equation 




0 


Equation (63) (along with eq. (54) in the interiors of the grids) allows 

(y n+ _ Q n ) Q k to assume its right value when the iteration process is carried to 

convergence (together with the corrective interpolation procedure described 
above). Both time accuracy and a spatial accuracy consistent with the order of 
interpolation used are maintained at the overlap boundaries with the use of equa- 
tion (63). 


Patch boundaries - The implicit patched-boundary condition (for the factored, 
implicit algorithm) is used along patch-boundary ABCD separating grids 3 and 4 
(fig. 54). The periodicity condition imposed on the upper and' lower boundaries of 
these two grids implies that the segments AB and CD need to be patched together in 
the calculation. While this leads to more involved programming logic, the procedure 
is conceptually straightforward. 
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Results 


In this subsection, we present results obtained for the rotor-stator configura- 
tion shown in figure 51. These results were obtained by integrating equation (54) 
and the boundary conditions described above. Four iterations were performed at each 
time step. The magnitude of the maximum residual in the system was reduced by more 
than an order of magnitude at the end of the fourth iteration. Approximately five 
cycles (a cycle corresponds to the motion of the upper boundary of grid 4 from its 
current position to the position occupied currently by the lower boundary of grid 4) 
were required to get rid of the initial transients and to establish a solution that 
was periodic in_time. The calculation was performed at a constant time-step value 
of .0.04 (the corresponding maximum CFL number in the grid system was approximately 
250). 

The dependent variables were nondimensionalized with respect to the ambient 
pressure (p^) and density (p m ). This yields 

u = M /y 

CO CO 

v =0 (inlet flow is axial) 

03 

where M m is a reference Mach number. For the purpose of initializing the depen- 
dent variables in the four grids and also specifying the Riemann invariants and 

R^j it was assumed that M m was equal to 0.1. However, because the quantities that 
are prescribed at the inlet boundary are R^, R^ and the inlet flow angle and not 
the dependent variables themselves, the average Mach number at the inlet (when the 
solution became periodic in time) was found to be approximately 0.08. The velocity 
of the rotor airfoil to was taken to be approximately 1.33 times the average inlet 
velocity (this is fairly close to the experimental value, which was 1.28). 

The Reynolds number for the calculation was chosen to be 100, 000/ in. The 
Baldwin-Lomax turbulence model (ref. 37) was used to determine the eddy viscosity. 
The kinematic viscosity was calculated using Sutherland's law. The viscous terms 
were evaluated only in grids 1 and 2. It was assumed that because of the sparseness 
of grid points in grids 3 and 4, the magnitudes of the viscous terms in these grids 
were negligible. 

In the figures that follow, there are several comparisons made with experimen- 
tal data. The following points must be kept in mind when evaluating these 
comparisons : 

1. The airfoil geometry used in the numerical calculation only approximates 
that used in the experiment (the actual configuration consisted of 22 stator air- 
foils and 28 rotor airfoils; the calculation was performed with 22 of each with an 
enlarged rotor geometry). 

2. The rescaling of the rotor geometry requires the modification of the Reyn- 
olds number to simulate equivalent conditions in the calculation, and it is not 
clear how this modification should be effected. 
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3. The inlet conditions were reasonably close to but not identical with those 
used in the experiment. Since the quantities that are specified at the inlet are 
not the dependent variables themselves, it is difficult to exactly specify experi- 
mental inlet conditions such as inlet velocity and inlet total pressure. The speci- 
fication of quantities such as inlet velocity may require an iterative procedure 
wherein the inlet values of the Riemann variables are continuously modified such 
that the time-averaged numerical inlet conditions such as velocity and pressure 
match the time-averaged experimental inlet conditions. 

4. The current analysis is only two-dimensional in nature whereas the actual 
flow is three-dimensional. 

5. The axial gap between the airfoils in the experiment was 15$ of the chord 
length. It is difficult to estimate the equivalent axial gap in the case of the 
modified rotor. The calculation was performed using an axial gap that was 15$ of 
the average chord length. 

Time-averaged pressures - Figure 55 shows the time-averaged pressure coefficient 
(Cp) as a function of the axial distance along the stator airfoil. The pressure 
coefficient is defined as 


C 

P 


p avg ~ (p t ) inlet 

(1/2)p. . .to 2 
inlet 


where P avK is the static pressure averaged over one cycle, (p^.) inlet i s ^he 
average total pressure at the inlet, and p^ n -j_ et is the average density at the 
inlet. Clearly, there is a good agreement between theory and experiment. A small 
separation bubble was found on the trailing edge circle of the stator in the numeri- 
cal results. This is seen as a sharp dip and rise of C toward the trailing edge 
on the pressure side. The experimental data also indicate such a separation. 
However, the magnitude of the pressure fluctuation obtained numerically may be 
suspect because the turbulence model was not tailored to yield accurate estimates of 
the eddy viscosity in such regions. 

Figure 56 shows the time-averaged Cp distribution for the rotor airfoil. The 
agreement is good except on the suction side of the rotor toward the trailing edge 
(4.0' < x < 7.0). This difference between experimental and numerical results is 
probably because of three-dimensional effects (caused by the low aspect ratio of the 
airfoil). A small separation bubble was found on the pressure side as in the case 
of the stator. The bubble is seen as a spatial fluctuation in pressure. Once 
again, the magnitudes of these fluctuations may be inaccurate. It is suspected that 
the increased mixing that exists in the bubble will result in larger eddy viscosity 
values and smaller pressure fluctuations in the real case (the sparsity of experi- 
mental points precludes any definite conclusions at the present time). 

Unsteady pressures - Figure 57 shows the magnitude of the pressure fluctuation 
C along the surface of the stator (plotted as a function of the axial distance). 
Tne quantity Cp is defined as 
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where P max and p m j_ n are the maximum and minimum pressures that occur over a cycle 
at a given point. The agreement between experiment and theory in figure 57 is 
fairly good considering that in the experiments the rotor moving past the stator is 
much larger than the actual one. Figure 58 shows the C distribution for the 
rotor. The comparison between experimental and numerica? values is not as good as 
in the case of the stator. However, the numerical calculation does predict all the 
qualitative features shown by the experiment. 


The periodic boundary conditions on the upper and lower boundaries of grids 3 
and 4 result in a periodicity in time of the solution within these grids. The 
pressure and other dependent variables become periodic in time after the initial 
transients dissipate. Any numerical scheme (including the interior and boundary 
schemes) that is devised to calculate rotor-stator flow must be able to simulate 
this periodicity in time. To test this particular capability of the interior inte- 
gration and boundary schemes, the pressure was monitored at selected locations on 
the stator. Figure 59(a) shows the locations at which the pressure was monitored, 
and figures 59(b)~59(d) show the variation of pressure over three cycles, as a 
function of time, at these locations. The periodicity in time is clearly observed, 
thus demonstrating the capability of the method in simulating periodic time- 
dependent flows. The magnitude of the temporal pressure fluctuations vary and was 
found to be maximum at the point B. A similar test was performed at selected loca- 
tions on the rotor. The test locations are shown in figure 60(a) and the temporal 
pressure fluctuations at these locations are shown in figures 60(b)-60(d). Clearly, 
the fluctuations are periodic in time. 


Figure 61 shows a three-dimensional plot of pressure as a function of time and 
position on the stator. To facilitate comparisons with a similar plot in refer- 
ence 30 it has been plotted with pressure decreasing in the upward direction (how- 
ever, the plot in ref. 30 does not contain the mean component of the pressure). 

This results in low-pressure regions appearing as peaks and high-pressure regions as 
valleys. The variation of pressure amplitude on the stator can be clearly seen in 
figure 61. The largest amplitudes occur on the suction side and close to the trail- 
ing edge. The amplitudes at the leading edge, and on most of the pressure side, are 
very small. There is a gradual increase in amplitude as we move from the leading 
edge to the trailing edge on the suction side. Another interesting feature is the 
sudden change in amplitude at the trailing edge as one moves from the suction to the 
pressure side. Figure 61 also shows the phase shift that occurs in the low pressure 
region L along the suction side. The low-pressure peaks close to the leading edge 
are almost 180” out of phase with respect to the peaks at the trailing edge. All 
the features mentioned above agree qualitatively with the experimental data of 
reference 30. 

Figure 62 shows the three-dimensional plot of unsteady pressures for the 
rotor. The pressures in this case are plotted with increasing values in the upward 
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direction. The largest amplitudes occur at the leading edge and on the suction 
side. The pressure amplitudes are rather small at the trailing edge and all along 
the pressure surface. The amplitudes decay as we move away from the leading edge on 
the suction surface. There is also an abrupt change in amplitude as we move from 

the suction to the pressure side at the leading edge. An interesting feature that 

can be observed in figure 62 is the phase shift that occurs in the high pressure 

region H as one moves away from the leading edge on the suction surface (the high- 

pressure peaks occur at different times). As in the previous case, the numerical 
data agree qualitatively with the experimental data. 

While it is difficult to make detailed comparisons of the unsteady pressure 
surfaces (figs. 61 and 62) obtained experimentally and numerically, it is a rela- 
tively simple matter to make comparisons of any one particular feature on these 
surfaces. Figure 63 shows the phase of the low-pressure region L (fig. 61) for the 
stator. The plot shows the time at which the low-pressure peak occurs as a function 
of arc length along the stator. The initial positions (t = 0) for the rotor and 
stator are shown in figure 51. The experimental values are below the numerical 
values on the pressure side and above the numerical values on the suction side. 

This is because the frequency at which the rotor airfoils pass the stator airfoils 
in the experiment is different than that in the calculation (there are fewer rotor 
airfoils in the calculation). This happens despite the fact that the tangential 
velocity- (w) of the rotor airfoils is the same in both the experiment and numerical 
calculation. While the numerical data can be modified to take into effect the 
differences in blade passing frequency (between the experiment and the calculation), 
truly accurate phase estimates can be obtained only with a multiple-airfoil 
calculation . 

Mach number contours - Figures 64-66 depict Mach number contours at various 
positions of the rotor with respect to the stator. Although the calculation was 
performed with only two airfoils, for the sake of clarity these contour plots depict 
several airfoils. The information regarding the additional airfoils is obtained 
from the periodicity condition. The thickening of the boundary layer as we move 
from the leading to the trailing edges of the rotor and stator is evident. The 
wakes associated with the airfoils are also clearly seen. In figure 64 the stator 
wake is below the rotor. In figures 65 and 66 an interaction between the stator 
wake and the rotor is observed. 

Figures 64-66 show contour lines that are smooth and continuous across the 
patch boundary. However, this is not the case with the contour lines in the vicin- 
ity of the overlap boundaries. It was found that the overlap boundaries gave rise 
to small-amplitude, high-frequency oscillations in the pressure and density values 
close to the overlap boundaries. The low Mach numbers in the calculation tended to 
sustain these oscillations because of the decoupling of the continuity, momentum, 
and energy equations. It is suspected that the overlap boundary procedure gives 
rise to these oscillations because conservation is not strictly enforced at these 
boundaries (the patch-boundary condition, on the other hand, is fully conservative). 
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Areas for Future Research 


There are several differences between the experimental conditions and the 
conditions under which the numerical calculation was performed. Chief among them is 
the difference in the number of rotor airfoils and the size of these airfoils. To 
accurately simulate the experimental data one needs to use a minimum of 1 1 stator 
airfoils and 14 rotor airfoils (in order to be able to implement the periodicity 
boundary condition at the proper locations). Such a calculation would be extremely 
time consuming. A better approach, and one that would justify the necessity for a 
multiple-airfoil calculation, would be to use 4 stator airfoils and 5 rotor air- 
foils; a total of 9 airfoils. The ratio of the number of airfoils (rotor vs. 
stator) for this configuration is 5/4 = 1.250, whereas the ratio in the real case 
is 28/22 = 1.273. Hence, the factor by which the rotor airfoils will have to be 
enlarged is 1.273/1.250 = 1.018 instead of the factor 1.273 used in the present 
calculation. The 9-airfoil calculation should yield better phase relationships and 
pressure amplitudes. 

The second problem that needs to be addressed is the turbulence model. Experi- 
ence indicates that an algebraic turbulence model does not accurately predict eddy 
viscosities in regions of separation and in wakes without requiring extensive "fine- 
tuning" of the model. Hence, it is important to use a more generally applicable 
turbulence model such as the two-equation models. Preliminary investigations show 
that the K - e turbulence model (ref. 38) yields more accurate estimates of the 
eddy viscosity in the separation regions near the trailing edges of the airfoils 
and, hence, results in smaller spatial pressure fluctuations. 

The conservative treatment of overlap boundaries is much more difficult than 
the conservative treatment of patch boundaries. For this reason, the overlap bound- 
aries in the present study were not made fully conservative. In the author's opin- 
ion this is the cause for the high-frequency oscillations in the dependent variables 
close to the overlap boundaries. This problem may worsen in situations where there 
are flow discontinuities passing through overlap boundaries. Hence it is preferable 
to either use only patched grids for multiple-grid calculations or to develop a 
simple way of treating overlap boundaries in a conservative manner. 

The improvements mentioned above are necessary in order to build a reliable 
rotor-stator analysis tool. Current research is therefore being focussed on build- 
ing an unsteady, Navier-Stokes, multiple-blade, patched-grid code which uses a two- 
equation turbulence model. 


SUMMARY 


To summarize, these lecture notes give a brief review of the work that has been 
done in the area of multiple grids, both patched and overlaid grids and then present 
a patch-boundary condition that was developed by the author. This patch-boundary 
condition is stable, accurate and easily applicable in generalized coordinate 
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systems. It is developed for both explicit and implicit integration schemes. In 
the case of implicit schemes, it is developed within the framework of an iterative 
implicit procedure that makes possible the decoupling of the calculations in the 
various patches (decoupling within a time step). The iterative implicit procedure 
is also included in the text. 

Results are presented for a variety of test cases. The example calculations 
include supersonic flow over a cylinder, blast wave diffraction by a ramp, the 
motion of vortices through patch boundaries, and subsonic flow over a cylinder 
demonstrating the feasibility of performing calculations on patches that move 
relative to each other. These calculations demonstrate the quality of solutions 
possible with the patch-boundary scheme in terms of both spatial and temporal 
accuracy. 

One of the most important applications of the patched-grid technology is solu- 
tion of problems that involve the motion of some parts of the system with respect to 
other parts. Typical examples are rotor-stator interaction in turbomachinery, 
helicopter rotor-fuselage interaction, propeller-nacelle interaction, and the coun- 
terrotating propfan. Two examples of rotor-stator interaction have also been 
included. The various boundary conditions employed to solve the rotor-stator prob- 
lem are discussed in some length. Results in the form of steady and unsteady 
pressures and phase plots are presented. Numerical results are compared with 
experimental results for one of the cases. A good agreement between experimental 
and numerical results is obtained. 

Current research is focussed on the following areas: (1) making the patch- 

boundary condition more accurate, and (2) extending the scheme to three-spatial 
dimensions. Preliminary results in both these areas of investigation have yielded 
encouraging results. 


48 



REFERENCES 


1. Rai, M. M. ; Chaussee, D. S.; and Rizk, Y. M. : Calculation of Viscous Super- 

sonic Flows Over Finned Bodies. AIAA Paper 83-1667, July 1983. 

2. Hessenius, K. A.; and Pulliam, T. H.: A Zonal Approach to Solution of the 

Euler Equations. AIAA Paper 82-0969, June 1982. 

3. Rai, M. M. ; Chakravarthy , S. R.; and Hessenius, K. A.: Metric-Discontinuous 

Zonal Grid Calculations Using the Osher Scheme. International Journal of 
Computers and Fluids, vol. 12, no. 3, 1984, pp. 161-175. 

4. Cambier, L.; Ghazzi, W. ; Veuillot, J. P.; et Viviand, H.: Une Approche par 

Domaines pour le Calcul d ' Ecoulements Compressibles . Cinquieme Colloque 
International sur les Methodes de Calcul Scientifique et Technique, INRIA, 
Versailles, France, 14-18 Decembre 1981. 

5. Beam, R. M. ; and Warming, R. F.: An Implicit Factored Scheme for the Compres- 

sible Navier-Stokes Equations. AIAA J., vol. 1 6 , no. 4, April 1978, 
pp. 393-402. 

6. Rai, M. M. : A Conservative Treatment of Zonal Boundaries for Euler Equation 

Calculations. AIAA Paper 84-0164, January 1984. 

7. Rai, M. M. : An Implicit Conservative Zonal Boundary Scheme for Euler Equation 

Calculations. AIAA Paper 85-0488, January 1985. 

8. Rai, M. M. ; and Chakravarthy, S. R.: An Implicit Form for the Osher Upwind 

Scheme. AIAA Paper 84-0088, January 1984. 

9. Rai, M. M. : A Relaxation Approach to Patched-Grid Calculations with the Euler 

Equations. AIAA Paper 85-0295, January 1985. 

10. Chakravarthy, S. R. : Implicit Upwind Schemes Without Approximate Factoriza- 

tion. AIAA Paper 84-0165, January 1984. 

11. Eriksson, L. E. ; and Rai, M. M. : A Stability Analysis of Various Patched Grid 

Interface Conditions for Hyperbolic Equations. To be published. 

12. Eriksson, L. E. : Euler Solutions on 0-0 Grids Around Wings Using Local Refine- 

ment. Proceedings of the sixth Gamm Conference on Numerical Methods in 
Fluid Mechanics, Gottingen, September 25-27, 1985. 

13. Usab, W. J.; and Murman, E. M. : Embedded Mesh Solutions of the Euler Equation 

Using a Multiple-Grid Method. AIAA Paper 83-1946-CP, July 1983. 


49 



14. Lombard, C. K.; and Venkatapathy , E.: Implicit Boundary Treatment for Joined 

and Disjoint Patched Mesh Systems. AIAA Paper 85-1503-CP, July 1985. 

15. Lombard, C. K.; Bardina, J.; Venkatapathy, E.; and Oliger, J.: Multi- 

Dimensional Formulation of CSCM — An Upwind Flux Difference Eigenvector Split 
Method for the Compressible Navier-Stokes Equations. AIAA Paper 83-1895, 
July 1983. 

16. Berger, M. J.: Adaptive Mesh Refinement for Hyperbolic Partial Differential 

Equations. Ph.D. Thesis, Department of Computer Sciences, Stanford Univer- 
sity, August 1982. 

17. Berger, M. J.: On Conservation at Grid Interfaces. ICASE Report Wo. 84-43, 

September 1984. 

18. Atta, E.: Component Adaptive Grid Interfacing. AIAA Paper 81-0382, January 

1981. 

19. Atta, E. H.; and Vadyak, J.: A Grid Interfacing Algorithm for Three- 

Dimensional Transonic Flows About Aircraft Configurations. AIAA 
Paper 82-1017, June 1982. 

20. Steger, J. L.; Dougherty, F. C.; and Benek, J. A.: A Chimera Grid Scheme. 

Mini-Symposium on Advances in Grid Generation, ASME Applied, Bioengineering 
and Fluids Engineering Conference, Houston, June 20-22, 1983. 

21. Lee, K. D.; Huang, M. ; Yu, W. J.; and Rubbert, P. E.: Grid Generation for 

Three-Dimensional Configurations. NASA CP-2166, October 1980. 

22. Yu, N. J.: Transonic Flow Simulation for Complex Configurations With Surface 

Fitted Grids. AIAA Paper 81-1258, 1981. 

23. Rai, M. M.: Navier-Stokes Simulations of Rotor-Stator Interaction Using 

Patched and Overlaid Grids. AIAA Paper 85- 151 9-CP , July 1985. 

24. Chakravarthy , S. R. ; and Osher, S.: High Resolution Applications of the Osher 

Upwind Scheme for the Euler Equations. AIAA Paper 83-1943, July 1983. 

25. Roe, P. L.: Approximate Riemann Solvers, Parameter Vectors, and Difference 

Schemes. J. Comp. Phys., vol. 43, 1981, pp. 357-372. 

26. Lyubimov, A. N. ; and Rusanov, V. V.: Gas Flows Past Blunt Bodies. NASA 

TT-F-715, February 1973- 

27. Shankar, V.; Kutler, P.; and Anderson, D. A.: Diffraction of a Shock Wave by a 

Compression Corner: Part II - Single Mach Reflection. AIAA J., vol. 16, 

no. 1, January 1978, pp. 4-5. 


50 



28. Srinivasan, G. R.; McCroskey, W. J.; and Kutler, P.: Numerical Simulation of 

the Interaction of a Vortex with Stationary Airfoil in Transonic Flow. AIAA 
Paper 84-0254, January 1984. 

29. Hessenius, K. A.; and Rai, M. M. : Applications of a Conservative Zonal Scheme 

to Transient and Geometrically-Complex Problems. AIAA Paper 84-1532, June 
1984. 

30. Dring, R. P.; Joslyn, H. D.; Hardin, L. W. ; and Wagner, J. H. : Turbine Rotor- 

Stator Interaction. J. Eng. for Power, vol. 104, October 1982. 

31. Benek, J. A.; Steger, J. L.; and Dougherty, F. C.: A Flexible Grid Embedding 

Technique with Application to the Euler Equations. AIAA Paper 83-1944, July 
1983. 

32. Pulliam, T. H., and Steger, J. L,.: On Implicit Finite-Difference Simulations 

of Three Dimensional Compressible Flow. AIAA J., vol. 18, no. 2, February 
1980, pp. 159-167. 

33. Steger, J. L.; and Sorenson, R. L.: Automatic Mesh-Point Clustering Near a 

Boundary in Grid Generation with Elliptic Partial Differential Equations. 

J. Comp. Phys., vol. 33, no. 3, Dec. 1979, pp. 405-410. 

34. Salas, M. D. ; Jameson, A.; and Melnik, R. E.: A Comparative Study of the 

Nonuniqueness Problem of the Potential Equation. AIAA Paper 83-1888, July 
1983. 

35. Pulliam, T. H.; and Steger, J. L. : Recent Improvements in Efficiency, Accu- 

racy, and Convergence for Implicit Approximate Factorization Algorithms. 

AIAA Paper 85-0360, January 1984. 

36. Rai, M. M. ; and Chaussee, D. S.: New Implicit Boundary Procedures: Theory and 

Application. AIAA J., vol. 22, no. 8, Aug. 1984, pp. 1094-1100. 

37. Baldwin, B. S. ; and Lomax, H.: Thin Layer Approximation and Algebraic Model 

for Separated Turbulent Flow. AIAA Paper 78-257, January 1978. 

38. Chien, Kuei-Yuan: Predictions of Channel and Boundary-Layer Flows with a Low- 

Reynolds-Number Turbulence Model. AIAA J., vol. 20, no. 1, January 1982. 


51 



TABLE 1.- COMPUTING TIMES FOR THE EXPLICIT AND IMPLICIT 

RELAXATION SCHEMES 


Type of scheme used 

Explicit 

Implicit 
relaxation (one 
iteration/step) 

Implicit 
relaxation (two 
iterations/step) 

Time per 100 steps, sec 

2.46 

3.00 

5.88 

Number of steps to 

2700 

280 

157 

converge 

Time to converge, sec 

66.42 

8.40 

9.23 


TABLE 2.- COMPUTING TIMES FOR THE FACTORED 
IMPLICIT SCHEME 


Number of 
iterations 

Computation cost parameters per step 


1 2 


Time per 100 steps, sec 4.19 8.24 

Number of steps to converge Unstable 94 
Time to converge, sec Unstable 7.75 
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"I- 4 


PATCH 

BOUNDARIES 



BOUNDARY ZONES BOUNDARY 


(a) An example of patched grids. 


Figure 1.- Types of multiple-grids used in finite-difference calculations. 




ZONE 2 ZONE 1 

(b) An example of overlaid grids. 
Figure 1.- Concluded. 
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Figure 2.- Zoning of the multiply connected region associated with three airfoils. 
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b METRIC DISCONTINUOUS GRID 




d DISCONTINUOUS GRID 


(a) Continuous. 

(b) Metric-discontinuous. 

(c) Metric-discontinuous. 

(d) Discontinuous. 

Figure 4.- Types of patched grids used in finite-difference calculations. 
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k = k max + 1 / 2 



k = k max> 


CM 


X 

ro 
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Figure 5.- Grid to illustrate global conservation property of conservative 

difference operators. 



ZONE 2 


THIS LINE 
CORRESPONDS 
TO m = 3/2, k = 1/2 


Figure 6.- Two-patch grid to illustrate patched-grid scheme in Cartesian 

coordinates. 
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Figure 



16.- Convergence history for the cylinder calculation (first-order 
explicit and implicit, pointwise relaxation calculations). 
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□ POINTWISE RELAXATION SCHEME 
(1 ITERATION/STEP) 

O POINTWISE RELAXATION SCHEME 



Figure 17.- Time-history of the maximum residual in the continuity equation 

(explicit and relaxation schemes). 
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A 


Figure 23.- Isobars after 24 integration steps (first-order factored, implicit 

scheme and five iterations/step). 
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CONVERGENCE HISTORY 



igure 24.- Convergence history for the cylinder calculation (first-order 
explicit and factored, implicit calculations). 
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EXPLICIT OSHER SCHEME 



Figure 25.- Time history of the maximum residual in the continuity equation 
(explicit and factored, implicit schemes). 
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ZONE 3 



Figure 30.- Grid for three-patch blast wave calculation. 
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ZONE 1 ZONAL ZONE 2 

BOUNDARY 


Figure 33.- Grid for two-patch vortex calculation. 



Figure 34.- Isobars for the vortex calculation (at initialization). 
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Figure 37.- Grid for single-patch vortex calculation. 



Figure 38.- Isobars for the vortex calculation (after 140 steps). 
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Figure 39-- Grid for two-patch cylinder (subsonic) calculation. 
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Figure 40.- Isobars at each step as the outer grid performs one rotation. 
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ZONAL BOUNDARY 

Figure 41.- Two-patch grid for the double-airfoil calculation. 
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Figure 44.- Pressure contours after 4.4 cycles. 
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Figure 46.- Pressure contours after 4.8 cycles. 
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Figure 47.- Pressure contours after 5.0 cycles. 
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STATOR BLADE 


INLET FLOW 
DIRECTION 

=^> ( 



DIRECTION OF MOTION 
FOR ROTOR 


SUCTION 

SURFACE 


Figure 50.- Rotor-stator geometry of Dring et al. (ref. 30). 




Figure 51.- Modified rotor-stator geometry. 
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STATOR BLADE 


ROTOR BLADE 



Figure 52.- "0" type grids for regions 1 and 2. 


SPACE FOR 
STATOR BLADE 



Figure 53.- Algebraically generated grids for regions 3 and 4. 
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OUTER BOUNDARY 


OF ZONE 3 



Figure 54.- Zoning of the rotor-stator problem showing patch and overlap 

boundaries . 



o PRESSURE SIDE, EXPERIMENT 
O SUCTION SIDE, EXPERIMENT 



LEADING TRAILING 

EDGE EDGE 


AXIAL DISTANCE, x 

Figure 55.- Time-averaged pressure distribution on the stator. 
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LEADING TRAILING 

EDGE EDGE 

AXIAL DISTANCE, x 

Figure 56.- Time-averaged pressure distribution on the rotor. 
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Figure 57.- Pressure-amplitude distribution on the stator. 
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O PRESSURE SIDE, EXPERIMENT 
O SUCTION SIDE, EXPERIMENT 
NUMERICAL RESULTS 



SUCTION LEADING PRESSURE 

SIDE EDGE SIDE 

< — AXIAL DISTANCE ALONG ROTOR SURFACE ► 

Figure 58.- Pressure-amplitude distribution on the rotor. 
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1.02r 


STATOR, POINT A 


(a) 

0 .5 1.0 

TIME, t 



(a) Stator locations at which unsteady pressures were monitored. 

(b) Unsteady pressure at point A on the stator. 

(c) Unsteady pressure at point B on the stator. 

'(d) Unsteady pressure at point C on the stator. 

Figure 59.- Unsteady pressures on the stator. 
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ROTOR, POINT A 



TIME, t 



. v/ - 

0 .5 1.0 0 .5 1.0 

TIME, t TIME, t 


(a) Rotor locations at which unsteady pressures were monitored. 

(b) Unsteady pressure at point A on the rotor. 

(c) Unsteady pressure at point B on the rotor. 

(d) Unsteady pressure at point C on the rotor. 

Figure 60.- Unsteady pressures on the rotor. 
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62.- Unsteady pressures 


on the 


rotor surface 



o PRESSURE SIDE, EXPERIMENT 
o SUCTION SIDE, EXPERIMENT 

NUMERICAL RESULTS 

PORTION OF NUMERICAL 

DATA DISPLACED BY ONE CYCLE 



Figure 63.- Phase of the low-pressure peak seen in figure 61. 



Figure 64.- Mach number contours (t = 0.00). 
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Figure 65.- Mach number contours (t = 0.25). 



Figure 66.- Mach number contours (t = 0.50). 
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